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FORTRAN IV PROGRAM FOR STUDYING SON - POLAR -MOLECULE COLLISIONS 

by R. Bruce Canright, Jr. , and John V. Dugan, Jr. 

Lewis Research Center 

SUMMARY 

Ion - polar -molecule collisions have been studied by numerically integrating their 
approximate equations of motion on a digital computer. In this report we present a 
FORTRAN IV program to do this integration (using conservation of energy as a step-size 
control) and to give results in the form of collision statistics, microfilm plots, and mo- 
tion pictures of trajectories. Included are the description of the ion -molecule system, 
the form of the input and output, the functions of the main program and subroutines, the 
COMMON structure, a sample case, and the complete FORTRAN IV listings. 


INTRODUCTION 
Problem Background 

The ion-molecule system has been studied extensively (refs. 1 to 3). By the use of 
a simple model, the Lagrangian can be obtained and the equations of motion derived (see 
ref. 4). Numerical integration with the aid of a computer gives instantaneous coordina- 
tes and velocities. The polar molecule can "look like" either a rod or a symmetric 
top. When the molecule is represented as a symmetric top, the laboratory coordinate 
system is as shown in figure 1 . 

The model potential energy has two terms, a permanent dipole contribution and an 
induced charge contribution. The Lagrangian for a symmetric-top molecule two scaled 
units long is written as 


L = T - V = — (X 2 + Y 2 + Z 2 ) + — (i 2 + r? 2 sin 2 £) + — (i// + r/ cos £) 2 

non 


+ 



(-X sin £ sin 77 + Y sin £ cos 77 + Z cos £) + 




where 


T model kinetic energy 

V model potential energy 

m reduced mass, MjMg/CM^ + Mg) 

X, Y, Z Cartesian coordinates of ion (see fig. 1) 

X, Y, Z Cartesian velocity components of ion 
| , V , xjy angular coordinates of dipole (see fig. 1) 

£ , ? 7 ,{[/ angular velocity components of dipole 

I^Ig moments of inertia (principal and about symmetry axis, respectively) 
p dipole moment of molecule 

e electronic charge 

r ion -molecule separation 

a electronic polarizability of molecule 



Figure 1. - Coordinate system for ion - polar-molecule pair. 
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Equations of Motion 


From this Lagrangian are derived 12 first-order coupled differential equations 
which must be integrated. These are the equations of motion. 


X s v„ = Jii- f-r 2 sin £ sin 77 - 3 XF ^ ^ae 2 X 


mr 


mr 


6 


Y *V 5 

mr 


= ( r 2 sin | cos 17 - 3 YF y j - 2ae 


mr 
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Z S t = j±. (r 2 cos £ - 3 ZF, 

z 5 V 
mr 


v)- 


2 oie 2 Z 


mr 


with F y = (-X sin £ sin 77 + Y sin £ cos 77 + Z cos £). 


X = V 


Y s V 


Z s V 


( 1 ) 

( 2 )- 

( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 


£ = (o t = 77 2 sin £ cos £ - (Z sin £ - X cos £ sin 77 + Y cos £ 

C O 






lo . o 

(^77 sin |+7? sin £ cos £) ( 7 ) 

! 1 


Me 


•0 = 0) = J=_ 

' 17 


(Y sin 77 + X cos 77) - 2 IjT7^ cos £ + ^ cos £) 

I| sin £ 


( 8 .) 


^ = t?| sin ^ - 77' cos £ 


( 9 ) 
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( 10 ) 


V = ^ (11) 

P s ( 12 ) 

Numerical Approach 

A theoretical check was made on the accuracy of integration: The instantaneous 
energy was compared at every time step to the initial energy, and the next time step was 
calculated based on this comparison. For this reason (a predictor-corrector check not 
required), a variable -order Runge-Kutta method was chosen (refs. 5 to 7), using fourth 
or fifth order for most problems. Details of choosing the (variable) time step for inte- 
gration appear in appendix B of reference 8. Essentially, the step size is inversely pro- 
portional to the relative difference between current total energy and initial total energy. 
This has proved to be much more efficient than simple step modification, for example, 
halving or doubling. 

Two important details of the numerical model are 

(1) The term sin(£) in the denominator of equation (8) goes to zero as the dipole 
lines up with the Z-axis, introducing a singularity into the equations of motion. This 
difficulty is handled by using two spherical (Euler when the dipole is a top, that is, \f/ ^ 0) 
coordinate systems for the dipole, an ,, unprimed M system and a "primed” system 

(ref. 1). The unprimed system (fig. 1) is used whenever the sin(|) term is greater than 
some chosen limit (the dipole is outside of some chosen polar cap). The primed system 
is used and the appropriate equations of motion are integrated whenever the dipole is 
within this cap. 

(2) The program provides an option of two collision potential terms: a hard -sphere 

potential V = °°, or a Lennard-Jones potential V ~ r - * 2 . Here r c is a chosen 
r r c 

collision distance (called CD in the program) . 


COMPUTER STUDIES 

Two kinds of studies can be made with this program: collection of statistics and de- 
tailed study of individual trajectories. For a collection of statistics , a large number of 
initial conditions must be generated "randomly. " Among other outputs, the program 
counts "captures; " by definition, the polar molecule captures the ion whenever the ion 
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comes within a certain prescribed distance from the molecule. An ion may be captured 
many times during the course of one trajectory; or it may be repulsed, never captured. 
The capture cross sections for various inputs can then be estimated (ref. 8). The pro- 
cedure is to integrate many sets of initial conditions for each of a number of values for 
the impact parameter b. Some fraction C R (b) of the number of sets for each b will be 
capture collisions. Then the capture collision cross section a is 

«c = ' jf° CrM d < b2 > 

where this integral is approximated numerically by using the calculated (b, C R (b)) points. 
Here the details of each trajectory are not important. 

For a detailed study of individual trajectories, the entire history of variables in an 
individual case is stored and plotted (refs. 9 and 10). Plots of relative translational mo- 
tion, rotational energy, and angles against the ion-dipole separation r, show the tra- 
jectory at a glance. Motion pictures of the system during its interaction provide a nat- 
ural view of the collision (in the center -of -mass system). (Plots and motion pictures 
were made on a Control Data Corporation Model 280 microfilm recorder; the arithmetic 
was done in double precision on an IBM 7040/7094 DCS.) The program presented herein 
is applicable for statistics, plots, and motion pictures. 


PROGRAM DETAILS 
input 


The input to the main program is of two types. First, there is the general input for 
a group of initial conditions. For convenience, this input is read in by means of 
NAMELIST (ref. 1), which includes these quantities: 


MASS 

MI 

MUE 

ALFE2 

HI 

ORDER 

TOTERR 


reduced mass of system, scaled units 
moment of inertia, I| 

Me, scaled units 
cte /2, scaled units 

initial separation, scaled to angstroms; placed in Q(l) 
order of Runge-Kutta integration, usually fourth or fifth 
total relative error in energy to be tolerated 
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maximum step -size increase factor 
maximum local error ratio to be tolerated 
initial step, scaled time 

sin(|), defines polar n cap, ” causes coordinate switch 

LOGICAL switch, .TRUE, if collision-time distance is input 

collision -time distance - the separation from which measurement of collision 
lifetime is started. (If not input, the program calculates 

CTD = 

where cr-^ is the Langevin cross section. ) Another criterion sometimes 
used (but not in the program as listed) is to start measuring collision time 
when interaction is ’’felt, ” that is, when the dipole rotational energy 
changes by 10 percent from its initial energy. The collision lifetime is the 
elapsed time the pair spends within this distance. 

collision distance (Every time the pair gets within this separation they col- 
lide, and a capture is counted.) 

maximum captures allowed for each trajectory 

boundary for Lennard -Jones polarization potential 

minimum step to be tolerated, scaled time 

^/Ip RM is zero for rod cases 

control array calling for still plots (If not all zero. ) 

■fVi 

program stores and plots every IPLOT point 

LOGICAL switch: if .TRUE. , truncates trajectories at 1000 points; if 
. FALSE. , uses disk storage 

Second, there is the set of initial coordinates Q and two labels, read in by 
(OP4F5. 3, 1P6E10. 3/OPF5. 3, 4X,2I3). Q(l) is not read in here, but is set to RI. 


Q(D 

r 

Q(5) 

1 

Q(9) 

V 

Q(2) 

e 

Q(6) 

R 

QUO) 

k 

Q(3) 

<p 

Q(7) 

e 

Q(ix) 


Q(4) 

V 

Q(8) 


Q(12) 



CD 

MAXCAP 

CD1 

HTS 

RM 

PLOT 

IPLOT 

CUTPLT 



MAXHR 

MAXER 

HI 

CHANGE 

CTDI 

CTD 
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Output 


For each case, the main program prints out the initial conditions, the final condi- 
tions, initial and final energies, impact parameter b, the initial velocity collision time 
and distance, various computing times, and the statistics of the step-size history (see 
also the section Sample Case). Units and scaling of the variables are given in table I. 


TABLE I. - UNITS AND SCALE FACTORS FOR VARIABLES 


IN MAIN PROGRAM 


Variable name 


Input units 


Output units 


Units inside program 


MASS 

MI 

MUE 

ALFE2 

HI 

CTD 

CD 

CD1 

HTS 

Q(l) or RI 
Q(2) 

Q(3) 

Q(4) 

Q(5) 

Q(6) 

Q(7) 

Q(8) 

Q(9) 

Q(io) 

Q(ll) 

Q(12) 

V 

B 


gxlO 20 

(g-cm 2 )xl0 36 

24 

(esu-cm)xlO 

(esu 2 -cm 3 )xl0 40 
1 4 

secxlO 

A 


gxlO 20 

(g-cm 2 )xl0 36 
(esu-cm)xlO 24 
|(esu 2 -cm 3 )xl0 40 ^ 


secxlO 

A 


14 


o 

A 

o 

A 

secxlO 14 

A 


rad 


o 

A 

o 

A 

secxlO 14 

cm 

rad 


gxlO 20 

(g-cm 2 )xl0 38 

0 A 

(esu-cm)xlO 

(esu 2 -cm 3 )xl0 40 

secxlO 14 

o 

A 

A 


O 

A 

14 

secxlO 1 

o 

A 


rad 


(cm/sec)xlO“ 8 

(rad/sec)xl0“ 14 


cm/sec 

rad/sec 


(cm/sec)xlO“ 6 

(rad/sec)xl0" 14 


rad 


t 

rad 

cm/sec 

cm 


rad 

cm/sec 

cm 


E or IE 
RE or IRE 
T 

REALT 

CT 

CTDS 

TRUN 

TPLOT 

CAPT 


ergxlO 8 

ergxlO 8 

sec 

sec 

sec 

cm 

sec 

sec 

rad 


ergxlO 

ergxlO 

secxlO 

sec 

sec 

cm 

sec 

sec 

rad 


8 

8 

14 
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Main Program 


The main program reads the input, sets up the case, monitors the integration, and 
prints the output. With only slight exceptions, it is the only routine which calls other 
routines; the call structure is very simple. The arithmetic and logic flow for this pro- 
gram are diagrammed in figure 2. 

The important switches local to this program are as follows: 


RS set to 2 if a capture occurs, otherwise set to 1 

CD1SWI set to 1 if CD1 > CD, signifying a Lennard -Jones potential (r - ^) term is to 

be included for this case; set to 3 if CD1 < CD 

CTSW keeps track of capture status 
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CD 



(b) Integration loop. 
Figure 2. - Continued. 



































(c) Output for a case. 
Figure 2. - Concluded. 


CD1SW initially CD1SWI; then keeps track, if necessary, of nearness to Lennard- 
Jones distance GDI 

REPEAT .TRUE, if bad step, . FALSE, if good 


Subroutines 

In addition to the main program, several subroutines are required. These sub- 
routines are summarized in table EL 






TABLE H. - SUBROUTINE SUMMARY 


Subroutine name 
(Entry) 

Purpose 

Called by- 

Calls on- 

JACK 

Evaluates the equations of motion at time t 

RK 


ENERGY 

Calculates various energy terms at time t 

Main 


RK 

Integrates from t to t + h, using 
variable -order, variable h Runge- 
Kutta method 

Main 

JACK 

SETRK 

Defines the Runge-Kutta coefficients for 
the current order 

Main 


PTUP 

Transforms dipole coordinates from 
primed to unprimed 

Main 


UPTP 

Transforms dipole coordinates from 
unprimed to primed 

Main 


SCT 

Transforms ion coordinates from spher- 
ical to Cartesian 

Main 


DARC 

Calculates a double precision cos -1 

Main 


^DS 

Makes plots; initializes 

Main 

HEAD; DOTS; BLANK, TEST 

DDR 

Stores points; entry in DDS 

Main 


DDP 

Plots points; entry in DDS 

Main 


HEAD or HEADM 

Makes headings 

DDP 


DOTS 

Draws dotted circles (for motion pic- 
tures) 

DDP 


BLANK, TEST 

Makes leaders (for motion pictures) 



PLTDUM 

Contains dummy entries for the purpose 
of this report (Entries are explained 
in comments cards, see section Pro- 
gram Listings.) 

Main; DDP 



a Two decks with same name; one used for plots, the other for motion pictures. 
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This COMMON area appears in the main program, and in DDS and HEAD; the variables 


are 


RE 

rotational energy 

V2 

velocity squared 

SXI 

sin(£) 

R 

ion -dipole separation 

T 

time 

KT 

rotational energy scale factor 

ID1,3D2 

plot labels 

Q(12) 

coordinate vector 

CLRP 

plot label switch 

RUN(3) 

plot label, BCD 

PLOT (5) 

plot selector array 

CUTPLT 

plot truncation switch 


(2) /JACKEN/EROC, RM,V, ETRC,EPLC, EPTC,XXC1,XXC2,XC, JCD1,EP,SIG, 
ETRAN, EPOT, EPOL: This COMMON area appears in the main program and in JACK 
and ENERGY; the variables are 

EROC Ij/2 

RM I 2 /I 1 

V velocity 

ETRC m/2 

EPLC ae 2 /2 

EPTC Me 

XXC1 Me/m 

XXC2 2ae 2 /m 

XC Me/I 1 

JCD1 switch for Lennard -Jones potential (If JCDl = 1, infinite barrier. ) 

EP Lennard -Jones interaction potential strength parameter 
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SIG 


characteristic Van der Walls diameter 


ETRAN translational energy 

EPOT potential energy 

EPOL polarization energy 




12 


where energy, E = EROT + ETRAN - EPOL - EPOT 

(3) /INTSET/C1, C2, C4, RN appears in the main program and in SETRK; the var- 
iables are a function of the order n of the Runge-Kutta method: 


Cl 

(2 n - l)" 1 

C2 

(2 n - 1) _1 

C4 

(MAXHR)' 

RN 

l/(n + 1) 


(4) /INTES/TOTERR, MAXER, HI, MAXHR, HMAX, HMIN, HAV, NS , NS F: This 
COMMON area appears in the main program and in SETRK; it contains the integration 
parameters 


TOT ERR 

MAXER 

HI 

MAXHR 

HMAX 

HMIN 

HAV 

NS 

NSF 


total relative error in energy to be tolerated 
maximum local error ratio to be tolerated 
starting time step for this case 
maximum step increase factor 
maximum step taken during this integration 
minimum step taken during this integration 
average step taken during this integration 
good step counter during this integration 
bad step counter during this integration 


(5) / RKC/NV , NVM 1 , C(11),B(12), A(66): This COMMON area appears in SETRK , 
and RK; these variables are the Runge-Kutta constants , constructed by SETRK and used 
by RK: 
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(6) /PCS C/ PCS appears in the main program and in DDS, JACK, and ENERGY. 

PCS = 1: not in polar cap, unprimed 

PCS = 2 in primed system 

(7) /LANG/ IS appears in the main program and in JACK and ENERGY. 

LS = . TRUE. : ju = 0, Langevin case 

LS = . FALSE. : not Langevin 

(8) /CIRCLE/RADIUS , RVIEW, ON appears in the main program and in DDS and DOTS 

(used for motion pictures only). 

RADIUS a function of collision distance and of distance from X-Y motion -picture 

plane; the radius of the dotted molecule or ion (Fig. 3 shows the function of 
RADIUS and RVIEW. ) 

RVIEW the scanning radius for the motion -picture frame, in angstroms (The ion and 
molecule are filmed when they are within two RVIEW of each other. ) 

ON If ON = . TRUE. , a dotted circle will be on the frame. If ON = . FALSE, no 

circle can be drawn. 


|£.2£7 

f.571 

J1 ... . 

§.87f| 
E. _ 


CS3- 


Hu»; 


too 20 

.9690. 6922. 480-4. 964E- 02-1 . 41 7E-04-1 . 937E-04-5. 54 IE- 02 8, 656E-Q2 l\ 




100 6 


>*£ 9E~ 03 £ • 99 1 E- 04- 1 . S36E- 04-1 . 917E-02-C . SS4E-C2 41 




.OT=l*CTBI = .F 


* *L3T= 1 j 1 » 1 > 1 j 1 jKT=. 259B- 1 * N9 :<CmP= 1 0$ 


*0NEj 


j 0 0 fl! 

■ in! 


10* C:-t:-r*1E=.=;* 001 = 0. DO, 00=9.00* HTS=l.E-6, RM=. OODO 

i i' i 1 h 1 r r ■ * * c 

.DC," 0RBER=4~TDTERR=. 0005D0”* ilftXHR=4 . D 0 » r«KER=27D07~ 

b ; I ■( » t . 8 

«ftSS^3T743^3 7hfl=T. 44Jp3» HUE=4.'9D^5 »PLFf 2*47490-37” 


54 5 5 S’| 
6|6 6 5 6| 
7(77 7 7 
3 8 8 8 8 
9 : 9 9 9 9 


FORTRAN STATEMENT I identification 

u‘ e s oo o j o o o r? o o o o o it p o u o o o o o o o o o o o c o odoToTFq o ooWbol FqTq oooooooaooa old ooCti FTo] 

6 I •• , : t; ij li • ! »<; ]. tl 11 ;. >> /i. . : ,a :s id jl JT » « l. A ]l ID ,1 «? « « « « 47 « 49 to SI */ SI » 55 M 57 M •» Kiel 62 6)M 63 bt 6f H .3 n tl n II H 13 It .. H J4 ao 

1 1 1 1 1 1 ! 1 tfi n u n 1 i 1 m 1 1 1 1 1 1 1 1 1 1 1 ! 1 1 f 1 m 1 1 1 1 1 1 1 m 1 1 1 1 1 1 11 1 1 1 1 nil 1 1 11 1 1 1 1 


2 2 2 2 2 2 2 2 7 2 2 2 2 ?? 2 2 ? 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 j 
3*3 3 3 3 3 1 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 I 
1 4 : 4 4 4 #|4| « <44444^! 4 44444 44 ^ 4444444444444444444444444444444444444444444444444444 44444| 
IS 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5-5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 
6|s 6666666566686668808688686666666666666666666666666686666666686666866886686 
V 7 7 1 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 ? 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 
|ajs 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 B 8 8 8 8 3 8 8 8 3 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 j 


23QE-01I 

; . ] 


Figure 3. - input card images. 
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Sample Case 


The input cards for a sample case are shown in figure 3. Two blank cards "trigger” 
the reading of the NAMELIST and then the RUN array only if plots are required; two 
Q- cards follow. This case produces the following output on the listing: 


$ONE 








MASS 

= 

0 • 2 74 CCCOCCGCCCOOCO— C2 » 

Ml 

0 . 1440000000 100000D-02, 

MUF 

C.4800CCCOOCOOOOC0D-O4 , 

ALFE2 

= 

0.449CCCCCCCCCCC00D-02, 

PI 

0.25 00000000 0000 00 C 02* 

ORDER = 

4, 

TQTERR 

* 

C . 50CQCCCCCCCCOCGCO-C3 » 

NAXHR » 

0.400c000000000000r 01, 

MAXER = 

C.2COOCCCCCCCOOOC09 01, 

HI 

= 

0.5C0CCCCCCCCCCCG0D Ol, 

CHANGE® 

5 • 0000000 E-Ol* 

CTCI = 

F , CTD 

= 0.1 rCCCCCOOOOOOOOPD 02, 

CC 

= 

0.2OOCCCOCCCCCCCCCO Cl* 

MAXCAP= 

10, 

CC1 

C. 

, HTS = 1.000iOnOE-06 

RM 

= 

0 . , 






PLCT 

= 

1. 

1 » 

1 , 

1, 


1, 

KT 

= 

Q.259CCC CCCOOOCCOQU-O 1 » 

IFLCT = 

1, 

CUTP LT= 

T, 


$ ENC 









CASE 1 

INPLT 

1.571 

2.2CS 6.073 

1.823 

-5.4290-02 

2 .99 ID— 04 

-i.esen-c* 

-1.91 

70 - C 2 -2.8840-02 

4.7?lD-0? 
5.870 100 



R 

THETA 

PHI 

ETA 

XI 

PS I 




INITIAL 

Q 

2. sod- c ? 

1 .5710 

2.2090 

6.0730 

1.8230 

5.870C 

V 

5. 5 CD 04 9 

4. 000— OQ 

CONG IT ION S 

CU 

-S.43D 04 

2.9SD 10 -1 

.86 0 10 

-1.920 12 

-2.880 12 

4.72C 12 

IE 

5.030-06 IRE 

8.470-07 

F INAL 

Q 

2. 51D-C7 

2.C671 

-0.7181 

-43.5264 

1 .5306 

26.4455 

E 

5.C30-C6 RE 

3.220-06 C 

VALUES 

OQ 

3.61D 04 

-2.78D 10 -1 

.910 10 

-6.20D 12 

-2.51D 12 

5 .450 12 

RATIT 

1.000 ti«e 

1.010-11 

COMP TIME 

e.87 

NS 

3 52 NSF 

67 

HAV 2 

.8 65 HM AX 

1C .44 £ 

EMIN 

C.C57 CAPT 

= 1.9065 

CCLL I S ION 

TIME 2 1 

. 12E-12 

DISTANCE 

8.23E-0 9 

RUN TIME 8.87 

PLOT TIME C . 




CASE 2 

INPLT 1.571 

3.965 C. 692 

2.480 

-4 .9 640-0 2 -1.417 C— 04 

- 1.5370-C4 

-5.54 10- 02 8.6560-02 

1.2800-01 









2.227 100 


R 

THE TA 

PHI 

ETA XI 

PS I 




INITIAL 

U 2.50D-C7 

1.5710 

3.9690 

0 .6920 2.4800 

2.2270 

V 

5.00D 04 B 

3 . 000— OR 

CONDITIONS 

CQ -4.S6U C4 

-1.420 10 -1 

. 94D 10 

-5.54D 12 8.660 12 

1.28D 13 

IE 

9. 64D-06 I RE 

6.2 3D— 06 

F INAL 

0 2.51D-C7 

C . 9857 

2.3194 

2.1625 0.2356 

1.5187 

E 

9.640-C6 RF 

7.360-06 C 

\ZfiLUES 

CQ 4.080 04 

1.430 C9 -2 

.78 0 10 

-3.74C 13 5.10D 12 

5.35D 13 

PATIO 

1.000 TIME 

9.400-12 

COMP TIME 

21.63 NS 

753 NSF 

256 

HAV 1.248 HM AX 

4.6 14 

hm in 

0.C59 CAPT 

= 2.3051 

COLL I S ION 

TIME 3.58E-12 

DISTANCE 

1.21E-07 

RUN TIME 21.63 

PLCT TIME 0. 
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In addition to this output, motion pictures or still plots can be made by including the de- 
sired subroutines. Figure 4 shows a sample motion -picture frame. Figure 5 shows 
sample still plots. 
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Rotational energy scaled Velocity scaled to initial 
to kT R » 0. 026 eV value of 5. 5x10^ cm sec'^ 



X axis 



(a-1) Relative velocity. 



Ion -molecule separation, A 
(a -2) Molecular rotational energy. 


(a) Variations of ion-molecule relative velocity and polar- 
molecule rotational energy during Ar + + CO single- 
reflection capture collision. 




Z axis 


<b) Variation of ion coordinates for NOj + NCI multiple- 
reflection capture collision. 


Figure 5. - Sample still plots. 
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(d-2) Azimuthal angle, (p ; vertical dashed line indi- 
dates passing through Oor 2tt. 


(d) Variations of coordinate angles for argon ion relative 
to CO molecule during capture collision with multiple 
reflections. 


Figures. -Concluded. 


Program Listings 

The complete program is listed below in FORTRAN IV (ref. 11). Not included are 
the routines which actually make the microfilm plots, because they are unique to this 
Center. A dummy deck, PLTDUM, explains these routines in comments cards so that 
others can be substituted. Of course, statistics were gathered, and even individual tra- 
jectories studied, without plots. 


19 




oonnonoooo n non noon n nonoooonoo on 


SIBFTC MAIN 


7 090. DOUBLE PRECISION SINGLE SYMMETRIC TOP I WITH ROD AS A SPECIAL 
C ASE I PROGRAM WITH 0080 PLOTTING OPTION 
'INTEGER PC SyRS.CTSW, ORDER, CASE »CLRP,PLGT,TNCOL, 

1 INCASE, CD IS W» CD ISWI 
LOGICAL LS 

LOGICAL REPEAT *C TO I »CUTPLT 

DOUBLE PRECISION TOTERR, MAXER *HI » MAXHR , HMAX, HMI N , H AV , 

1 C 1 ,C2 yC4 yRNy 

2 EROC ,RM, V, ETRC »E PLC «EPT C* XXC1 ,XXC2*XC«CDi, 

3 RE y V2 y SXI *Ry TyKT.Qy 

4 R I , MASS , MI yMUE y ALFE2 , 

5 B y CAPTyCQ yRDyCTDyEylEyIREy RAT I Oy RAT IQL» 

6 H,EA,EL,ER,TMAX,HN,Ql,SZ2yTl,Z«PSID 
DOUBLE PRECISION X ♦ SORT, ARCOSy EP .SI G . ETRANy EPOT , EPOL 
LOGICAL ON 

EXTERNAL JACK 
DIMENSION Z( 12) ,01(12) 

COMMON FOR MOVIE DECKS CNOT STILL PLOTS) 

COMMON/CIRCLE/ RADI US.RVIEWyON 

COMMON FOR INTEGRATION PARAMETERS AND STATISTICS 

COMMON /INTES/ TOTERR , MAXERyHL, MAXHR .HMAX.HMIN, HAV . NS. NSF 

COMMON FOR COMMUNICATION WITH SETRK. QUANTITIES ARE FUNCTIONS OF T 
INTEGRATION FORMULA ORDER 
COMMON /l NT SET/ C1,C2,C4,RN 

COMMON FOR TRANSMITTING CONSTANTS TO ENERGY AND JACK (JACK « DIFFE 
ENTIAl EQUATION EVALUATOR) THE TRANSLATIONAL VELOCITY, V, ALSO 
APPEARS. 

COHMON/JACKEN/ EROC y RM, V, ETRC ,.EPLC ,EP TC , XXC1 . XXC2, XC, J C01 , 

1 EP , SI G. ETRAN, EPOT.EPOL 

COMMON FOR L ANGEVIN SWITCH. LS*.TRUE. IFF MUE*0. 

COMMON/LANG/ LS 

COMMON HOLOS PRIME D- UNPRI MEO COORDINATE SYSTEM SWITCH 

1 = UPRIMED 

2 » PRIMED 
COMMON /PC SC/ PCS 

COMMON FOR TRANSMITTING CONTROL PARAMETERS, HEADING INFORMATION, 
AND PLOT QUANTITIES TO PLOTTING SUBROUTINE 

COMMON /DDC / REyV2ySXI,RyTyKTylDlyID2.Q(12) . CLRP. RUN! 3 ) . PLOT ( 5 ) , CJ 
IT PLT 

DATA C A S£ /0/,L ORDE R/0/ , REPEAT/. FALSE. /,CTD/L.D1/, CTO I/. FALSE./ 

DATA NCOL . TNCOL .NCASE ,TNCASE/0 ,0,0,0/ 

THESE VARIABLES ARE READ IN VIA GENERAL INPUT (NAMELIST) CARDS. 
MASS, MI, MUE, AND ALFE2 MUST BE SPECIFIED. IF CTDI IS SPECIFIED 
AS TRUE, THEN CTO MUST ALSO BE SPECIFIED. IF THE ENTIRE PLOT 
VECTOR IS NOT SET TO 0, THEN RUN AND KT MUST BE SPECIFIED. 

MASS REDUCED MASS OF ION-DIPOLE PAIR 
MI MOMENT OF INERTIA IN SCALED UNITS 

MUE DIPOLE MOMENT FACTOR IN SCALED UNITS 

ALEE2 POLARIZATION FACTOR IN SCALED UNITS 
R I INITIAL SEPARATION IN ANGSTROMS 
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ORDER ORDER OF RUNGE-KUTTA FORMULA TO BE USED 

TOTERR TOTAL ERROR TO BE TOLERATED OVER INTEGRATION 

MAXHR MAXIMUM STEP SIZE INCREASE FACTOR 

MAXER MAXIMUM ERROR RATIO TOLERATED 

HI INITIAL STEP SIZE IN SCALED TIME UNITS 

CHANGE S IN ( XI )=CHANGE INITIATES COORDINATE TRANSFORMATION 

CTO I LOGICAL VARIABLE - COLLISION TIME DISTANCE INPUT SWITCH 

CTD COLLISION TIME DISTANCE IREAD IN ONLY IF CTDI = .TRUE- ) 

CD REFLECTION DISTANCE (CAPTURE DISTANCE FOR C.S. CALCS) 

MAXCAP MAXIMUM NUMBER OF CAPTURES PERMITTED 

CDI BOUNOARY OF POLARIZATION POTENTIAL, USED AS L— J SWITCH 

HTS MINIMUM STEP SIZE ALLOWED (IN SCALED TIME UNITS) 

RM RATIO OF SECOND MOMENT OF INERTIA TO FIRST MOMENT OF INERT 

FOR ROD CASES RM=O.DO. 

PLOT PLOT CONTROL VECTOR, IF(PLOT{ I ) .NE.O) THE I TH PLOT WILL B 
DRAWN 

KT AVERAGE ROTATIONAL ENERGY 

I PLOT PLOT EVERY IPLOT TH POINT 

CUTPLT LOGICAL VARIABLE - IF TRUE, PLOTS ONLY 124-9 POINTS TO AVOI 
TAPE DELAY 

NAMELIST/ONE/ MASS,M I ,MUE,ALFE2, RI» ORDER, TOTERR, MAXHR, MAXER, H 
1 1, CHANGE, CTDI, CTD, CD, MAXCAP, CDI, HTS, RM, PLOT, KT, I PLOT, CUTPLT 

REDEFINE FUNCTIONS AS DOUBLE PRECISION 
SQRT(X)=DSQRTIX) 

ARCOS( X)=DARCOS( X) 

PRESET GENERAL INPUT QUANTITIES 
R I = 25, DO 

EP, SIS L ENNARD-JONE S (CUSP ) POTENTIAL CONSTANTS 

EP=6.88D-6 
S I G=3« 500 

0RDER=4 

TOTERR = .00100 
MAXHR = 4, DO 
MAXER = 2. DO 
HI = 2, DO 
CHANGE 3 . 707 
CD = 1 . DO 
MAXCAP = 10 
HTS= 1, E-6 
RM=Q.DO 
CTOI=. FALSE. 

I PLOT= 1 
CUTPLT=. TRUE, 

DO i 11=1, NPLTS 
PLOT ( 1 1 ) *0 

SKIP TO TOP OF PAGE 
WRITE (6,31) 

READ A CARD 

READ ( 5,32) (QII), 1=2,12), lDl,ID2 

R DOT NONZERO INDICATES A CASE SET 
IF (Qf 6). NE.O. ) GO TO 4 
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THETA ZERO INDICATES BLANK SET - NEXT SET IS GENERAL INPUT 
IF (Q( 2J.EQ.0. I GO TO 3 

OTHERWISE' COMPUTE AND PRINT OUT COLLISION PERCENTAGES*, (SKIP IF NO 
CASES HAVE BEEN RUN YET. ) 

IF INCASE. EQ, 05 GO TO 2 

PCC0L = 100.*FL0AT(NCDL ) /FLOAT INCASE) 

TNCOL=TNCOL+NCOL 

TNCASE=TNCASE+NCA$E 

TPCC0L=100. ‘FLOAT I TN COL ) /FLOAT ( TNC AS E ) 

WRITE (6,33) NCOL.NCASE, PCCOL, TNCOL, TNCASE , TPCCOL 

NCOL=0 

NCASE=0 

GO TO 2 

READ GENERAL INPUT SET AND WRITE OUT ALL GENERAL INPUT VARIABLES 
READ (5, ONE) 

WRITE I 6, ONE) 

PASS RADIUS OF CORE TO MOVIE DECKS (IF ANY) 

RADIUS = CD 

READ IN 18 CHARACTER RUN LABEL IF THERE IS TO BE ANY PLOTTING 
IF I (PLOT! 1) .NE.O) .OR. I PLOT! 2) .NE.OJ.OR. I PLOT I 3) .NE.O).OR. (PLOT! A) 
1 .NE.O) .OR* (PL0T15) .NE.O) ) READ (5,34) RUN 

COMPUTE CONSTANTS FOR ENERGY AND DIFFERENTIAL EQUATION EVALUATION 

SUBROUTINES 

ETRC = • 5D0*MASS 

EPLC = . 5D0»ALFE2 

XXC2 = 2.D0*ALFE2/MASS 

SPECIAL CASE IF MUE = 0 
IF I MUE. NE.O.DO ) GO TO 50 
LS = .TRUE. 

EROC = 0. DO 
EPTC = 0.00 
XC = D.DO 
XXC1 = 0. DO 
RM = O.DO 
GO TO 51 

50 EROC = . 5D0*MI 
LS = .FALSE. 

EPTC=MUE 

XC=MUE/MI 

XXC1=MUE/MASS 

51 CONTINUE 


IF RUNGE-KUTTA FORMULA ORDER HAS BEEN CHANGED, CALL SETRK TO 

CALCULATE NEW COEFFICIENTS 

IF (ORDER. EQ.LORDER) GO TO 2 

LORDER=ORDER 

CALL SETRK (ORDER) 

SET CUSP POTENTIAL SWITCH INITIALIZATION VARIABLE 

m i <:ui a -s 

IF ICD1.GT.CD) CDISW 1=1 
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RETURNS TO READ A NEW SET OF CARDS 
GO TO 2 

NEW CASE - COUNT IT AND CLOCK IN 
CASE=C ASE + 1 
CALL TIMLFT (TINS 
NCASE=NCASE+1 

WRITE OUT THE CASE CARD 

WRITE (6*35) CASE, (Q( I) , 1=2, 12 > , IDl, ID2 

INITIALIZING CALL TO PLOT ROUTINE 
CALL DOS 

SET CUSP POTENTIAL SWITCHES 
CD I SW=CD 1 SW I 

EP.SI3 LENNARD— JONE S l CUS P ) SWITCH SETTING 

JCD1=I 

RETRIEVE INITIAL CONDITIONS AND UNSCALE THEM FOR OUTPUT 
Q( 1 )=R I 

Q(6) = -DABS ( Q ( 6 ) ) 

DO 5 1=2,6 

z c i )=a( i ) 

5 Z< 1+5) = Ql I+5)*1.D14 
Z(I) = Q(l)*l.D-8 
Z ( 6 ) = Q(6)*1.D6 
Z( 12) = 3( 12 ) 

COMPUTE VELOCITY, IMPACT PARAMETER, ESTIMATED INTERACTION TIME, AN 
COLLISION TIME DISTANCE 
SZ2 = DS I N ( Z ( 2 ) ) 

V=SORT( Z(6)*Z(6)+Z(1 ) *Z(1)*(Z{7)*Z(7)+(Z(8)*SZ2)**2)) 

TMAX = 2.D6+RI/V 

B = (Z( 1 ) **2*SQRT ( Z( 7)**2+(Z(8)*SZ2**2)**2) )/V 

IF l .NOT .CTDI ) CTD = SQRT( ALFE2/ (ETRC*V*V) )*l.D-2/B 

OUTPUT INITIAL CONDITIONS, VELOCITY, AND IMPACT PARAMETER 
WRITE (6,36) (Z(I) ,I=1,5),Z(12),V,B,(Z(I ),I=6,ll> 

TRANSFORM TRANSLATIONAL VARIABLES TO RECTANGULAR COORDINATES 
CALL SCT (01 1)»Q(6)»0(1) ,Q(6)»1> 

STORE INITIAL CONDITION VECTOR 
DO 6 1=1,12 

zm=o( i ) 

INITIALIZE TIME, STEP SIZE MAXIMUM AND MINIMUM STEP SIZE, 

REPULSION SWITCH, STEP COUNTER, BAD STEP COUNTER, AND CAPTURE 

COUNTER 

T = 0» DO 

HM I N=H I 

HMA X=H I 

HN=H I 

RS= 1 

NS = 0 

NSF = 0 

NCAP=0 
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PERMUTE INTERGRATION VARIABLES AND DERIVATIVES TO ORDER R, THETA»P 
ETA ? XI ? PS I 
PSID=3S 11 ) 

00 7 K = 6, 10 

I= 16 -K 

Q(I+1!=Q(I) 

QS6)=Q( 12) 

Q( 12 )=PSID 

SET PRIMED-UNPRIMED COORDINATE SWITCH TO UNPRIMED 
PCS=i 

COMPUTE INITIAL ENERGY 
CALL ENERGY! I E , IRE »Q 11) » QI 7 ) ) 

INITIALIZE ENERGY RATIO AND R 
RATIOL = 1 * DO 
R=R I 
RETIRE 

RECORD INITIAL VALUES FOR PLOTTING 
CALL DDR 

OUTPUT INITIAL ENERGY AND ROTATIONAL ENERGY 
WRITE (6*37) IE, IRE 

INITIALIZE COLLISION TIME SWITCH AND COLLSION TIME 

CTSW=1 

CT=0» 


BEGIN A STEP - SET STEP SIZE AND TRIAL VARIABLE 

H=HN 

Tl= T 


BRANCH IF WE ARE IN PRIMED COORDINATE SYSTEM 
IF ! PCS.EQ.2 ) GO TO 9 

WE ARE IN UNPRIMED SYSTEM. BRANCH IF WE SHOULD STAY IN THIS SYSTEM 
IF (ABStSINI SNGL ( Q ( 5 ) ) )). GE„ CHANGE > GO TO 10 

TRANSFORM TO PRIMED COORDINATE SYSTEM 
PCS=2 

CALL UPTP (Q(I)*Q(7),0m«Q(7)) 

BRANCH TO INTEGRATION CALL 
GO TO 10 

WE ARE IN PRIMED SYSTEM. CHECK TO SEE IF WE SHOULD STAY IN PRIMED 
SYSTEM. BRANCH TO INTEGRATION CALL If WE SHOULD 
9 XI = ARC0S(-DSIN(Q(5 ) )*DSIN(Q{4) ) ) 

IF (ABSISINlXl ) )„LT. CHANGE) GO TO IQ 

TRANSFORM TO UNPRIMED COORDINATE SYSTEM 
PCS = I 

CALL PTUP (Q(l) ,QI7) ,Q(1),Q!7) ) 

CALL RUNGE-KUTTA INTEGRATION SUBROUTINE. ORDER OF FORMULA HAS 
ALREADY BEEN SPECIFIED® 

10 CALL RK(12,T1,H,Q,Q1, JACK) 
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CALL ENERGY ROUTINE AND COMPUTE ENERGY RATIO 
11 CALL ENERGY! E»RE»Q1( l ) «Q1( 7) ) 

RAT 1 0=E / I E 

COMPUTE ERROR MADE IN THIS STEP 
EL = DABS l RAT IO-RATI OL ) 

COMPUTE ERROR ALLOWED FOR THIS STEP 
EA — H* ( TOTERR - DAB S { RATI OL-l . DO ) I/TMAX 

COMPUTE RATIO OF ACTUAL ERROR TO ALLOWED ERROR 
ER=EL/EA 

TURN STEP REPEAT SWITCH ON IF ERROR RATIO IS TOO LARGE 
IF (ER.GT.MAXER) REPEAT= .TRUE. 

BRANCH TO MAXUMUM STEP SIZE INCREASE SECTION IF ERROR IS VERY SMAL 
IF (ER.LT.C4) GO TO 21 

COMPUTE SIZE OF NEXT INTEGRATION STEP 
HN = H*(1.D0/ER)**RN 

ENTER HERE WITH NEW STEP SIZE 

2 CONTINUE 

TRY TO PREVENT HANGUP1ESP. HCL ) 

I F ( EL. GT. I .D-6 .OR. .NOT .REPEAT) GO TO 60 
REPEAT = .FALSE. 

60 CONTINUE 

BRANCH TO STEP FAILURE SECTION IF REPEAT SWITCH IS ON 
IF (REPEAT) GO TO 22 

SUCCESSFUL STEP - COUNT IT, RECORD LAST RATIO, ADVANCE TIME, 

UPDATE DEPENDENT VARIABLE VECTOR, UPDATE MAXIMUM AND MINIMUM 
STEP SIZES, RECORD LAST R, COMPUTE NEW R. 

NS=NS+ 1 
RAT IOL=RATIO 
T=T 1 

DO 13 1=1,12 

3 0(11*3111) 

HMAX=DMAX1(HMAX»H) 

HMIN=DMIN1(HMIN»H) 

R=SQRT (Q{ 1 )»Q( I)+Q(2 ) *Q ( 2)+Q(3)*Q( 3) ) 

FORK ON COLLISION TIME SWITCH 
GO TO I 14, 15,16) ,CTSW 

CHECK TO SEE IF WE HAVE PASSED COLLISION DISTANCE ON THE WAY IN. 

IF SO, RECORD TIME AND RESET SWITCH TO SEARCH FOR COLLISION DISTAN 
ON THE WAY OUT. 

4 IF (R.GT.CTD) GO TO 16 
CTSW=2 

CT 1 =T 
GO 10 16 

CHECK TO Sfct IF WE HAVE PASSED COLLISION DISTANCE ON THE WAY OUT. 
IF SO, RECORD TIME, COMPUTE ELAPSED TIME IN SECONDS, AND RESET SWI 
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C TO BYPASS THIS SECTION. 

15 IF ( R» LT» CTD ) GO TO 16 
CTSW=3 
CT2 = T 

CTMCT2-CT1 )*1.E-14 
C 

C RECORD EVERY I PLOT TH POINT FOR PLOTTING 

16 I F ( MOD ( NS » I PLOT ) . EQ* 0 ) CALL DDR 
C 

C CUSP POTENTIAL SWITCH SETTING SECTION 

GO TO 117,18,19) ,CD1SW 

17 IF (R.GT.CD1 ) GO TO 20 
CD1SW=2 

JCD 1 =2 
GO TO 19 

18 IF {R.LE.COl ) GO TO 19 
CD 1 SW= 1 

JCD1=1 
GO TO 20 
C 

C BRANCH IF CAPTURE CONDITION HAS OCCURRED. IA CAPTURE HAS OCCURRED 

C IF R IS LESS THAN THE COLLISION DISTANCE AND IS DECREASING, 1 

19 RD= ( Q( 1 )»Q(7)+0I2)*018)+Q(3)*Q<9) >/R 

IF ( (R. LT.CD ) .AND. (RO.LT.O.DO) ) GO TO 27 
C 

C BRANCH TO REPULSION SECTION IF CASE IS FINISHED. (A CASE IS FINISH 

C IF R IS GREATER THAN THE INITIAL R=RI) 

20 IF (R.5T.RI) GO TO 24 
C 

C RETURN FOR THE NEXT STEP 

GO TO 8 
C 

C MAXIMUM STEP SIZE INCREASE SECTION 

21 HN=MAXHR*H 
GO TO 12 

C 

C STEP FAILURE SECTION - TURN OFF REPEAT SWITCH, COUNT THE FAILURE, 

C AND RETURN TO REPEAT THE STEP UNLESS NEW STEP SIZE IS TO SMALL. 

22 REPEAT=. FALSE. 

23 NSF=NSF+1 

IF ( ABS ( SNGLI HN) ) .GE.HTS) GO TO 8 
C 

C STEP SIZE TOO SMALL - WRITE OUT SIGNAL, REDUCE CASE COUNT, SET PLO 

C CASE TYPE LABEL SWITCH, AND BRANCH TO OUTPUT SECTION 

RATERR = ABS(SNGLIRATIOL-l.DO) ) 

WRITE 16,38) RATERR 
NCASE-NCASE-1 
CLRP=4 
GO TO 28 
C 

C FORK ON REPULSION SWITCH 

C 1 = REPUSI ON 

C 2 = CAPTURE HAS ALREADY OCCURRED 

24 GO TO «25,26),RS 
C 

C REPUSI ON SECTION - WRITE OUT SIGNALS® SET PLOT CASE TYPE LABEL SWI 

C AND BRANCH TO OUTPUT SECTION 
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25 WRITE (6,39) 

CLRP=2 

GO TO 28 
C 

C COUNT CASE AS CAPTURE, WRITE OUT NUMBER OF CAPTURES THIS CASE* 

C SET PLOT CASE TYPE LABEL SWITCH, AND BRANCH TO OUTPUT SECTION, 

26 WRITE (6,40) NCAP 
NC0L=NC0L+1 
CLRP=I 

GO TO 28 
C 

C CAPTURE SECTION - COUNT CAPTURE, BRANCH TO CAPTURE CASE END SECTION 

C IF MAXIMUN NUMBER OF CAPTURES PERMITTED HAS OCCURED. 

C OTHERWISE TRANSFORM TRANSLATIONAL VARIABLES TO SPHERICAL 

C COORDINATES, CHANGE SIGN OF R DOT ( RADIAL VELOCITY), TRANSFORM TRAN 

C LATIONAL VARIABLES BACK TO RECTANGULAR COORDINATES, SET REPULSION 

C SWITCH TO INDICATE CAPTURE, AND BRANCH BACK TO RESUME INTEGRATION, 

27 NCAP=NCAP+1 

IF (NCAP.EQ.MAXCAP) GO TO 26 

CALL SCT (0( I),Q(7),Q(1) ,Q(7),2) 

Qm=-am 

CALL SCT (Q(1),Q(7),Q(1) ,Q(7)»i) 

RS=2 
GO TO 8 
C 
C 

C OUTPUT SECTION - CASE IS FINISHED. 

28 CONTINUE 
C 

C IF ROTATIONAL VARIABLES ARE IN PRIMED COORDINATE SYSTEM, TRANSFORM 

C TO UNPRIMED SYSTEM. 

IF (PCS.EQ.2) CALL PTUP (Q ( 1 ) , Q( 7 ) , Q ( I ) , Q ( 7 )) 

C COMPUTE SCATTERING ANGLE 

CAPT=ARCOS ( (Q(7)*Z(7)+Q( 8) #Z ( 8 > +Q ( 9 ) *Z ( 9 ) ) /SORT ( ( Q (7) *Q ( 7 ) +Q ( 8 ) *Q ( 
I8)+Q(9)*Q<9) )*(Zm*Z(7)+Z(8)*Z(8)+Z(9>*Z(9>n ) 

c 

C TRANSFORM TRANSLATIONAL VARIABLES TO SPHERICAL COORDINATES. 

CALL SCT (Q(1),Q(7),Q(1) ,Q(7),2) 

C 

C UNLESS THERE HAVE BEEN NO SUCCESSFUL STEPS, COMPUTE AVERAGE STEP 

C SIZE AND STORE IT FOR INITIAL STEP SIZE NEXT CASE 

IF (NS.EQ.O) GO TO 29 
HAV = T/DBLE(FLOAT(NS) ) 

H I =HAV 
C 

C COMPUTE RATIO OF FINAL ENERGY TO INITIAL ENERGY 

29 RAT 1 0= E/ I E 
C 

C UNSCALE INDEPENDENT AND DEPENDENT VARIABLES 

T = T* 1 . D- 14 
Q( 1 ) = Q(l)*l.D-8 
Q t 7 ) = 0(7) *1.06 
DO 30 1=8, 12 
30 Oil) = Oil ) * 1 ® Dl 4 
C 

C READ CLOCK AFTER INTEGRATION AND BEFORE PLOTTING 

CALL TIMLFT ( TBFP ) 
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c 

C PLOT THIS CASE 

CALL DDP 
C 

C CLOCK OUT AND COMPUTE COMPUTER TIME FOR THIS CASE 

CALL T I ML FT ITOUT) 

REALT=(T IN-TOUT )/60* 

C 

C WRITE OUT FINAL CONDITIONS? ENERGY? ROTATIONAL ENERGY* ELAPSED 

C COMPUTER TIME? AND INTEGRATION STATISTICS* 

WRITE 16*41) IQ<I)?I = l?6)?E?RE?IQm? 1 = 7? 12) ? RAT 1 0 ? T? REALT ? M$*MSF? 
1HAY?HMAX?HMIN 
C 

C WRITE OUT SCATTERING ANGLE 

WRITE i 6? 42 ) CAPT 
C 

C SET COLLISION TIME TO CASE TIME IF COLL I SON DISTANCE IS NOT LESS 

C THAN INITIAL SEPARATION 

IF (CTD a GE*RI) CT=T 
C 

C SCALE COLLISION TIME DISTANCE 

CTDS=CTD*KE-8 
C 

C COMPUTE RUN TIME AND PLOT TIME 

TRUN=CTIN~TBFP)/60« 

TP LOT= ( T8FP-TOUT )/60« 

C 

C WRITE OUT COLLISION TIME? COLLISION TIME DISTANCE? INTEGRATION TIM 

C AND PLOTTING TIME 

WRITE { 6? 43 ) CT?CTDS?TRUN?TPLOT 
C 

C RETURN FOR ANOTHER CASE 

GO TO 2 
C 
C 

31 FORMAT ilHl) 

32 FORMAT i 0P4F5, 3 ? 1P6E 10 . 3 /0PF5* 3 , 4X* 21 3 ) 

33 FORMAT C///6H CASES? I6?11H COLLISIONS? I6?8H PERCENT? F6« 1 / 1 2H TOTAL 
1 CASES? 16? 17H TOTAL COLL IS IONS? 16? 8H PERCENT ?F6®1///1 

34 FORMAT I3A6) 

35 FORMAT I////5H CASE? I 5 ? 3X? 5H INPUT ? 5X ? 0P4F7* 3 ? 1P6D12 * 3/ U2X ? 0PF7* 3? 
1 2X? 2 1 5 ) 

36 FORMAT I24X? 1HR? 8X ?5HTHE TA , 7X? 3HPHI * 8X? 3HETA * 8X? 2HX I ? 7X? 3HPS I / IX ? 7 
1HINITIAL?7X?1HQ?3X?1PD11*2® 0P5FII® 4?4X? 1HV • 2X • 1PD1 1 . 2?4X ? 1H8 ? 2X? IP 
2011*2/ IX? 10HC0ND ITI0NS?3X? 2HDQ? 3X? 1P6D11-2) 

37 FORMAT « 1H+, 87X?2HIE « 2X, lPD11.2»2X«3HIRE f 2X« 1PD11.2) 

38 FORMAT f 1 2 7X? 1HS * 1PE9 * 2 ) 

39 FORMAT I I23X ? 9HREPUL S ION ) 

40 FORMAT 1 123X « 1HC * 16) 

41 FORMAT { 1H+?5HFINAL? 9X?IHQ?3X?1PDI1*2?0P5F11«4?4X?1HE?2X?1PD11*2?3 

1X?2HRE?2X? 1PD11*2/1X? 6HVALUES?7X?2HDQ? 3X?1P6D11«2?4X?5HRATIO?OPF9« 
23? IX?4HTIME?2X*IPD11 «2/lX?9HC0MP T I ME? 3X ? 0PF7 * 2 ?4X* 2HNS? 3X ? I 5? 4X ? 3 
3HNSF?3X?I 5 ?4X? 3HHAV? 3X?F8® 3 ? 4X ? 4HHMAX?3X?F8® 3? 4X?4HHMI N? 3X? F8® 3) 

42 FORMAT I 1H+? 106X?6HCAPT =?F8*4I 

43 FORMAT C IX ? 14HC0LL IS I ON TI ME ? 1PE 12* 2 , 6X? 8HDI ST ANCE , IPE 12 * 2 » 6X? 8HRU 

IN TIME?0PF7*2? 6X? 9HP LOT TIM£?0PF7«2) 

END 
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SUBKUU FINE JACK ( T , X,D X) 

DiMENS IUM X{ IB) h)X( 12 ) 

IN T EGER P c s 
LOGICAL L 3 

OUUBLt PricC 1 SION I ,X,DX,EROC ,RM,V,ETKC,EPLC,EPTC,XXC1,XXC2,XC, 

1 VLJ, R,R2 »K3 ,R5,R6 ,X1 ,ETA,SXi »CET A»SETA, 

2 CXI , P 31 » UP S l sOETAfDXI , VPOT , XXC 3 , E P , S I G 

3 , 1 l'K AN , LPUT, EPOL 
COM MU N /PC SC/ PCS 
COMMON /LANG/ LS 

COMMON /JACKEN/ EROC ,RM , V , t T kC » E PLC ,EPTC , XXC i , XXC2 , XC, JCD1 , EP» S IG 
1 , E TRAN, E POT, E PUL 
DO 1 1 = 1,6 
i Dxm = xi i+o) 

ETA = X<4) 

XI = XI5) 

PSI = XI o) 

UXI = DXI 3) 

L)E f A = 0X14) 

UP SI = OXIo) 

R = 0 jwRTIXI 1)*XI1) + XI2)*XI2) + XI3)*X(3J) 

k 2 = R *R 

R 3 = R *R 2 

R3 = R2*R3 

SETA = OSINIETA) 

CETA = DCUgIETA) 

SXI = OSIN(XI) 

CXI = GC0S1XI) 

GO lu (20,30) » PCS 

20 VPOT = 3.00*1 X( I ) * SXI * SB T A - XI2 ) *3X1 *C£T A * X(3)*CXI) 

OX I 7) = XXo 1* (R2* 3XI*SETA - Xll)*VP0T)/R5 
DX( 8) = -XXC1*(R2*3XI*CETA + X( 2 )*VP 0 n/R 5 
0 X 19 ) = X Xc 1 * ( R 2*C XI - X(3) * V PO T) / R 5 

24 R6 = R 2*R 3 

XXC 3 = XXC2/R6 
00 23 1=1,3 

25 UX( 1 + 6 ) = 0X( IH1-XU ) *XXC3 
GO TU I 2o , 2 7 ) , JCU1' 

27 VLJ=EP*( SIG/R )**12 
DO 29 1=1,3 

29 DX( 1 + 6 ) = 0X( I+o) + 1 2 . D 0* X { I )*VLJ/R2 
2o IHLSi GO TO 40 

UXI 10) = ( XC* ( X( 1 )*CETA + X(2)*SETAI/R3 - 2 . 00 *DET A*UX I*CX I * 

1 RM*DX 1*1 DP3I+DE TA*CXI ) ) /SXI 

DX ( 11) = 0ETA**2*SXl*CXI + XC * ( C XI * ( X ( 1 ) *SET A-X( 2) *CET A) - X( 3 )* 

1 SXU/R3 - RM*SX1*I'UPSI+DETA*CXI )*UETA 
DXI 12 ) = DETA*D XI * 6X1 - DX(10)*CXI 
GO TO 50 

30 VPOT = 3. DO* ( X ( 1 ) *C XI - X (2 ) * SXI *CE TA - X ( 3 1 *SX I *S ET A I 
DXI 7) = XXC 1*1 R2*C XI - XI l) *VP0T)/R5 

DXI 8) = -XXC 1*IR2*GXI*CETA 4- X ( 2 ) * VPOT ) / R5 
DXI 9) = -XXC1*(R2*SXI*SETA + X I 3 1 * VP GT ) / R5 

34 R6 = R 3*R 3 
XXC 3 = XXC2/RO 
DO 33 1=1,3 

35 UXI 1 + 6) = UXIi+6) - XII ) * X XC 3 
GO TO ( 36 , 37 ) , JCD 1 
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3? VLJ=EP*f SIG/R)**12 
DO 39 1 = 1(3 

39 OX I 1 + 6 i = DX4 1+6) 4- 12.00*X4I UVLJ/R2 
36 IF! LSI GO TO 40 

OXiiOJ = IXC*4X(2)*SETA-X(3)*CETA)/R3 - 2. DO*D£T A*DXi *CX l * 

1 RM*D X I* 4 DP S t+DE TA*C X l ) ) / SXI 

DXCI11 * D£TA**2*SXI*CXI - XC*( X41) *SXI+CXI*IX(2)*CET A+X43 )*SETAI) 
1 /R 3 - RM*SXI*(DPSU-OETA*CXI ) *DE TA 
0X4 121 = DETA*D XI* SXI - 0X410)*CXI 
GO TO 50 

40 DX( 10) - O.DO 
DX4 in * O.DO 
DXi 12) = O.DO 

50 RETURN 
END 


$IBFTC ENRGY. 

SUBROUTINE ENERGY4E.ER0T, X,DX) 

INTEGER PCS. CLRP, PLOT, RUN 
DOUBLE PRECISION E.EKOT.X.DX, 

1 EROC »RM»V»ETRC»EPLC»EPTC»XXC1»XXC2»XC» 

2 RE,V2,SXI,R»T»KT»Q,k2,GXI»EP»SIG» ETRAN, 

3 XI .ETA.PSI »DXI .DETA.OPSI ,EPOL , E POT , CD12 
DOUBLE PRECISION VLJ 

LOGICAL CUTPLT 
DIMENSION X i 6) «DX(6) 

COM MON /PC SC/ PCS 

COMMOM/DDC/ RE.V2.SXI .R.T.KT .101.102. 0412) ,CLRP,RUN(3) . 

* PLOT! 5), CUTPLT 

COMMON /JACKEN/ EROC ,RM, V, £ TRC ,EPLC ,EPTC ,XXC1 , XXC2 , XC, JCD1 » EP. S I G 
1 .ETRAN. EPOT ,EPOL 
ETA = XU) 

XI - X4 5) 

PSI = X4 6 1 
DETA = DXI 4) 

DXI = DX4 5) 

DPSI = 0X46) 

SXI = DSINI XI) 

CXI = DC0S4 XI) 

R = DS0RT4 X( 1)*X(1) + X42 )*X(2) «• X43)*X<3)) 

R 2 = R*R 

EROT = EROC* UOXI*DXI * 4DETA*SXI )**2) + RM*f OPSI+OETA+CX l )**2) 

V 2 » 4 0X( 1 )*DX( 1) + OX(2)*DX(2 ) DX43)*DX(3) ) 

GO TO 41,2). PCS 

1 EPOT = ( EPTC /( R*R2 ) ) * ( X 11) *SXI *DSl N< ET A) -X 4 2 ) *SXI *DCOS { ET A ) + 

1 XI 3)*CXI ) 

GO TO (4,5), JCD1 

2 FPOI = <EPTC/IR*R2) )*(X(l)*CXl - X42) *SXI*DCOS I ETA) - XI3) *SXI* 

1 DSINIETAI) 

GO TO 44, 51, JC01 

4 EPOL = EP LCy <R2*R2 ) 

GO TO 6 

5 VLJ=EP*4SIG/R)**12 
£POL=EPLC/4R2*R2)-VLJ 
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6 t Its rt.Sl = t- V* i. # V f 

F =t-Ki i r + 1 I K AN -tH.il. -t.pt |T 
K f r UK M 
t NU 


$ IBFTC RK . 

SOBROU I 1 'I E KK. I Nt » TjH i X r XF ,F ) 

DOUBLE PRECIS 1 UN T »H , X ,XF ,A , a ,C »G, Y , 0 TE MP 
DIMENSION XINE),XF(NE) 

DIMENSION LAI 11 ) ,GI Id ,12) , Y( 181 
COMMON /KKC/N V,NVMl ,C 1 ID ,B 112) ,Alo6) 

DATA LA /0,1 ,3 ,6,10, 15 ,21 ,28 ,36 ,43 ,55/ 

CALL FIT, X,0< l»i>) 

DU 30 L = i , NVM 1 
DU 20 1 = 1 » NE 
DTEMP = 0 « DO 
M=L A( L ) 

DU 1U K=t,L 
M =M + 1 

10 DTtMP =OTEMP+A IM)*G (I ,M 
20 Y I I ) =X ( I ) *H*0 TEMP 

CALL FIT + H*C (L) » Y »G 1 1 »L+l) ) 

30 CONTINUE 

OU SO 1 = 1, NE 
DTEMP = O.UO 
DO 40 L = 1 , N V 

40 DTEMP =0 TEMP +B(L)*G( l ,L) 

30 XFi I)=X( 1 )+H*DTEMP 
T = T+H 
RETURN 
END 


$ IBFTC SETRK . 


SUBRGU I IN E SETRKINMI 

DOUbLE PRECISION A,B,C »C1 »C2 »C4 , RN, 

1 TOTERRfMAXEKjHl , MA XH R » H M A X » H M I N » H AV * 

2 TS1,CN,CD,AN,AD,BN,BD 
INTEGER NTH 8) ,NT2 ( 8),NT3( 7) ,NTI B),NVA( 7) 

INTEGER P 

OIMEN S ION CN( 36) ,CD (3 6) , AN (14 8) , AD (3 6) ,BN(43 ) , BO (7 ) 

COMMON /RKC/ N V »N VM 1 »C 111) ,8(12) , A I 66 ) 

CUR MON /I N TSE T / Cl ,C2 ,C4 , KN 

COMMON /INTES/ TO TtRK s MAXER , HI , MAXHR ,HMAX »HM1 N ,HAV , NS » NSF 
DATA NVA 72,3,4,6,7 ,9,12/ , 

1 NT 70,1 ,2, 3 ,4, 3 ,6,7 / , 
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2 NT 1 /i,2,4» 7,12,18 ,26,37/, 

3 N T2 /l, 3, 6, 10, L6, 23, 32, 44/, 

^ N T3 /1,2,5,1L,26,47,83/ 

DATA CN 

1/ 1.00, 1.0 0, 1. DO, i. 0 0, 1.00,1. 00,1.00 ,2.00 ,1 .00 ,2. DO, 4. DO, l. DO, 2. DO, 

21.00, 1 «D 0 , 1« j 0 j 1.00,2.00,1.00,1. GO 1 1 « DO , 8 * DO , 1 . 00 , 5 » 00 , 1 » 0 0 , 1 * DO , 

31.00, 1.00, 1.0 0, 1.00,1.00,2. 00, 1.00,5. 00,5. DO, 1.00/ 

DATA CO 

1/ l. DO, 2. JO, 1.00,2.00,2.00 ,1.00,3 .00,5,00 , l .DO ,3.00 ,5.00,3. 00,3 .DO, 
23*00, 2.00,2.0 0, i.O 0,9. 00, 3. 00 ,2.00 ,6.00 ,9 .00 ,9.00,6.00, 1.00, 9.00, 

36.00. 4.00. 10.00.6.00.2.00 .3. 00.3.00 .6 .00 .6.00 . 1.00, 

DATA BN 

1/ l. 00, 1.00, 1.00, t. 00, 1.00,1.00,2.00 ,2. UO , l .00 ,23.00,0.00, 125.00, 

A 0.00,-81 .00, 125.00, 

211.00, C.JC, 8 1.0 0,8 1.0 0,-32. 0 0, -32. 00, 11. DO, 110 2 0 1.00, 0.00, O.DO, 

376 79 3b. 00, 63 5 04 0.0 0,-5 9049. 0 0, -5 9049. 00, 63 50 A 0.00, 1 10 20 1 . DO , 41 . 00, 

40.00, 0.00, o. JO, 0.00,216.00,2 72.00,27.00 ,27 .00 ,36.00, 180.00, 41. 00/ 
DATA 80 

1 /2.U 0,6. DO, 6. 00, 192. 00, 12 0.00, 2 140320. DO, 840. 00/ 

DATA AN / l.UU, l. DO, -1. 00,2.00, 1.00,0. DO ,1.00 ,0.00 , O.DO, l. DO, l. 00, 

1 4.00, 6.00, 1.00,- 12.00,15.00,6.00,90.00,-50.00 , 8 . 00 ,6 . 00, 3 6 .00 , 

2 1.01, 6.00,0.00,1.00,0.00,2.00,1.00,4.00,-1.00,-1.00, 18.03,-3.00, 

3 -6.00,0.00, 9.00,-3.00,-6.00,4.00,9. DO ,-36 .03 ,63.00 ,72 .DO, 0.00, 

4 -64.0 0, 2.0 0, 1.00,3.00,1.00,0.00 ,3.00 ,23.00 ,3.00,21.00,-8. DO, 
5-4136.00, 0.00,-13 584.0 0,5264.00,13104.00,10513 1.00,0.00,33 2016.00, 

6 -107 744.00,-2842 56.0 0, 17 01 .00,-775229. DO ,0. 00 ,-2770950 .03 , 

7 17 35 l 36. U 0,25472 16. 3 0,8 1891. 00, 328536. 00 ,23869.03,0.00,-122304.00 

8 ,-20364.00,695520.00,-99873.00,-466560.00,241920. 00, l.DO, 1.00, 

8 3.00, 1.0 0,0. DO, 3. DO, 2 9. 0 0, 0.00, 33. DO, -12. 00, 33. 00,0.00,0. 00,4.00, 

9 125. 00,- 21. 00, O.DO, O.DO, 76. 00, 12 5. DO ,-162. 03 , -30 . 00 , 0 . 00, 0 . DO , 

A - 32.0 0, 12 5.00,0.00, 99.00,11 75. 00,0.00 ,0 . 00 , -345b. 00, -6253 .00, 

B 8424.00, 242.00,-27.0 0 ,2 93.00 ,0. 00 , 0 . DO , -852 . DO ,-l 3 75 . 00, 1 836 .00, 

C -1 18.00, 162. 00, 324. 3 0,1 3 03. 0 0,0. Du, 0.00 ,-42 60 .00, -68 75. DO, 999. Dl, 

0 1030.00, C. 00, 0.0 0,162. 00, -85 95. 00, 0.00, O.DO, 30 720. 00, 4875 O.DO, 

£ -66086.00, 376.00,-729.00,— 1 944 .00,-1296.00,3240.00/ 

DATA AO 

1 /l. 00, 2. DO, 1 .00, 2.00,2.00,1.00,3.00,25.00,4.00,81.00, 75. DO, 3. 00, 

2 3.00, 12. 00, 16. 00, 8. JO, 44. 00, 9. 00, 12. 00 , 8. 00, 216. 00 , 729. 03 , 

3151632 .00 . 1375920.00.251888.00 .9.00 .24.00 . 16.03 .5. 02. 972.00. 3, .00, 
4 243.00,324.00, 324.00 , 162.01 ,4428.00/ 

N N = N M 

I F ( NN .£0.0) Kfc TORN 

IF( MAX HR. L£. O.DO) MAXHR = 4.00 

P =NN+ L 

TSi - i.OO#*NN 
C 1 = 1.00/(781-1. DO) 

C 2= T S 1*C 1 

C4 = 1 .UO/(MAXHK**P) 

RN = 1 .00/0BL£(FL0AT(P) ) 

N =N T ( N N ) 

N V=NV A (N ) 

NVM 1=N V- 1 
I 1=NT 2 (N ) 

I2=NT2(N4 n- 1 
K =0 

Du 10 1=11,12 
K =K 1 
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10 b(K) = bM ( 11 / B D C N ) 
1 1=NT UN) 

I2=NT UN* 11- I 
K =0 
L =0 

00 20 1 = 11,12 

K =K * 1 

C(K 1 = CNUI/UU) 
DO 20 J=L,K 
M=N T3( N l + L 
L =L *1 

20 A(l } = AN(M)/AD< 1 1 
K fc r UK N 
END 


SIBFTC P TUP . 

SUbkUU UME P TUP ItP ,L»t P ,E ,Dt ) 

DUUBL £ PRECISION E P ,0E P , E , DE , 

1 y£ ,OUE »yEP»UQEP ,P,T,6,PD»TD,SD, 

2 PP,TP,SP,PPU,rPU,SPO, 

3 SNP,CP»ST»Cr»CrP,STP,SN» 

4 SPP ,CPP ,SSP,CSP 

DOUBLE PRECISION SI N ,COS, AKCC US , ATAN2 , X , Y 

0 HEN SION yE( 6) ,0yE(6) »OEP ( o) ,D QE P 1 6 1 , E (6 I , DE < 6 ) , E P ( 6 } , DEP < 6 I 
EQUIVALENCE <yb (4) »P) , <yE (6) , T ) , (QE (6) ,S) , (DyE (4J , PCI , ( DUE ( 5 I , TO), 
KDQEI 6!,SD),(WEP(4J ,PP) , (yEP(5) , TP 1 , <yEP(6) ,SP) , < OyEPi 4 I , P PD ) , 
2(Dy£P ( 51, TPD) , { DQ£P(6l,bPD) 

S IN t X) = OSIN(X) 

CO Si X ) = OCOSl X) 

AKCCOSUJ = DAkCuS(X) 

ATAN2 ( X, Y ) = UATAN2(X,Y) 

DU 18 1=4,6 
UEP ( I I =EP ( I I 
18 DUEP ( I )=DEP( I ) 

SPP =S INI P P I 
CPP =COS(PP I 
STP =S1N( TP I 
CTP =CU Si TP I 
S SP =S I N( SP ) 

CSP =CQS( SP 1 
C T=- S TP# SPP 
T=ARCCUS(CT) 

st=sin( n 

SN = OSIGNl i.uo,sn 
P=ATAN2(C TP*SN, STP#CPP*SN) 

SNP =S INi P ) 

CP = CUS(P I 

S=A TAN 2( l -CSP*CPP+C IP* SPP* S SP ) * SN, ( S SP*CP P+CT P*S PP*CS P I *SN 1 

SO=U SPD* C TP + PPO ) * SNP- ( rPO*iPP-SPO*SI P*CPP) *CP )/ST 

TD = i TPG*SPP-SPD*STP*CPP*SD*$T*CP!/6NP 

PD=-( TPO* CPP* SPO* STP* SPP+ SD*C T I 

UO 22 1=4,6 

E ( 1 1 =JE( 1 1 


33 



22 0 E ( I j =Uwt ( I ) 
00 30 1=1 f 3 
Ei I I =E Pi 1 i 
30 DEi I ) = 0EP { i) 
RETURN 
END 


$I6FTC UPTP« 

SUBROUTINE UP TP ( E f OE » fc P » OE P I 
DOUBLE PRECISION E*OE ? EP 9 DtP 9 

1 QE y OUt , QE P ?D wE P » P » T f S » PO , T 0 , S.D , 

2 PP f TP tSP fPPDjTP0 9 SP0f 

3 SNP , C P ? S T fCT ? S S ?C S ♦ C T P , S T PtSNt 

4 SPP,CPP 

DOUBLE PkECISION S l N ,C US t ARCC OS * AT AN2 » X * Y 

DIMENSION Qt( 6) ,JgE (6) ,QEP( 6) ,DQEP(6) ,E(6) f DEI 6) , EPI6 ) , OEPI 6 ) 
EQUIVALENCE (QEi*) ,P) *{QE<5) , T) , IQE(6i P SI * i OQE I 4 ) , PD) f ( DQE(5),TD)* 
Li OQci 6 1,301, l QE Pi 4 I f PPi * (QEP(5 J , TP J , ( QE P i 6 US P J f I OQEP ( 4 ) * P PO J f 
2I0QEP { 5) , TPO) , (OuEPibl , SPDi 
SlN(X) = DSINiX) 

cosi x ) = ocosi x) 

AKCCUS(X) = OAKCUSiXi 
AT AN 2 ( X, Y ) = UATAN2(X,Y) 

DO lb 1=4,6 
QEi H=E« I J 
lb OQE { i)=OE( U 
3NP = bliN(P J 
CP=C0 SIP ) 

ST=SIN( n 
CT =C0 S ( T ) 

SS= 31 N (SI 
CS=CUS(S1 
C TP =ST*SNP 
rp=ARCCQb(crp i 
STP=SIN( TP ) 

SN = OSi^NI WOO, STPi 
PP=ATAN2(-CT* SN, ST*CP*SN) 

SPP = 3 I NiPP ) 

GPP =COS(PP ) 

SP = AT AN2( (CS*uP-CT*3NP*S:>) * SN i SS*CP*CT *S NP*CS ) *S N ) 

SP 0 =i ( $0* $ T*CP- TO* SNP ) *CPP-( SD*CT+ PD) *SPP)/STP 
TPD = { TD*SNP-SO*ST*CP+ SPD* $ T P&C PP i / SPP 
PPO = r D*CP *SO*ST*SNP-SPO*CTP 
DO 22 1=4, b 

EP ( I J =QEP ( I 1 
22 DEP i I) =UQEP( i ) 

00 30 1=1,3 
EP ( I ) = E( I ) 

30 DEPU) =DE( I) 

KETUkN 

ENO 
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$ IBFTC SC T 


SUBROU 11NE xnS.DbtL »DC tMuDE ) 

uiMtN siun sm »dgi n ,c ( n »oc i n 

UGUBL E PRECISION S?UG?u?DC?K?DRjSr?CT?DTR$SP?DPR? 
l X*QXfY$uYyZjDZ?K2?CP 

GO TO 110,20) s MODE 
1C R = SI l) 

OK = 0 Si 1 ) 

ST = OSIN ( S( 2 ) ) 

CT = OCO S ( SI 2) ) 

DTK = K*JSI2) 

SP = DS1NISI3)) 

CP ■= OCOSl S( 3 I i 
DPR = k*U SI 3) 

Cl U = R* ST*CP 
C( 2 ) = k* SI*SP 
Ci 3 ) = k*CT 

DC ( U = OR*gT*CP * 0 TR*C T*CP - D PR^S T *SP 
UCIZ) = OR*Sr*SP + DTR*CT*SP + DPR*bT*CP 
OCI 3) = DK*C T - J TK*SJ 

GO TO 30 
2C X = Ci 1) 

dx = oci l ) 

Y = C( 2) 

DY = DCI2) 

L = CIS) 

OZ = OCI 3 ) 

R 2 = X*X ♦ Y*Y + L*L 
SI 1 ) = DSUR T I K 2 ) 

Si 2) = OARCUSIZ/Bl 1) ) 

SI 3) = UA TAN2 I Y , X ) 

OSI 1) = I X*DX + Y*0 Y + Z*DZ)/S<1) 

OSI 2) = I Z*J3(l)-S<l)*oZ)/{S(l)*USwRTtR2-Z*Z) ) 

DSI 3 ) = I X*UY-Y*DX)/IR2-Z*Z) 

30 RETURN 
END 


$ IBFTC DARC. 

DOUBLE PRECISION FUNCTION DARCUSIX) 

DOUBLE PRECISION X, Y , P 1 »H A LF P I * DA TA N» OSulRT 

OAT A PI/3. I Al 3926335 8979320 0/ »H ALF P I / l . 37 079o 3 26 794897 00/ 

IF! UABSIX ) .GT. 1 .DO) CALL ARERRI3 6HOARCOS CALLED WITH ARGUMENT 

1 l.i) 

ifi x. Eg. o.oo) go ru 10 

Y -= UA IAN IDSORT 11.00- X^X) / X) 

IF l Y .L I.O.DO) Y = Y «- PI 
DARCOS = Y 
GO TU 20 

10 DARCOS = HAL F P 1 
20 RE I URN 
END 



u-i r 


$ IBE TC DOS » * 


SUtJRU U T IN L DUS 

LOGICAL NuPLU I ? Sfv I H E T 9 C U T P L T 

INTEGER pcs.clrp.plut 

DOUBLE PRECISION 0 , V2 ,K »KE » T , KT , SXX 

COMMON /OJo / kE*V 2* SXX*R»T»Kr f ID1 ,[D2 ,Qtl2) ,CLRP, 

1 RUN< i),PLOT(5) *CUTPLT 
COMMON /PC SC/ PCS 

REAL RAJ 1000) , V A ( 1000) »REA( 1000) S XXA( 1000 ) * Y Y A ( 1000 ) , Z Z Al 100U ) » 

1 THETAAJ 1000) » XAJ 1000) , YAt 1000) »ZA (1000 ) » GAMMAAJ 1000 ) , PHI AJ 1000 ), 

2 0 A( 1000 » 12) , 1 V< 2) »K TC (2 ) »NEGP 

EQUIVALENCE 1 oA J 1 , 1 i , RA J 1 ) ) » <QA(1 ,2) ,VAlll ) , ( QA ( 1 , 3 ) » REA ( 1 ) ) t (Q 

1AJ 1,4) »X*Ai 1) ) , <Qa< 1 >5) ,YYA(1) ) » (QA(1»6) ,ZZA(l) ), ( QA( l, 7), THETA 
2A( 1) ) , IJAI 1, 8 ) tXAJ 1) ) 9 (QA (1 ,9) ,YAll) ) , ( wA ( 1 , 10 ) , Z A ( 1 ) ) * (QA(1,1 
31) » GAMMAA ( 1) ) , ( QA ( 1 » 1 2 ) , PH I A { 1 ) ) 

DATA NPA/12/,LREC /I 000/ ,R TO/5 7. 2 95 7 8 / , EOP/K/ , NEOP/O./ , 6EG/4H$R IB/ 
l» EN0/4H$X IE/* 1V/12H /»KTC/12H- /» LRECM 1/999 / 

DATA NPL r S/5/ 

NUPLO T = • TRUE . 

00 1 11=1 » NPL T S 

1 if- (PLOT! 1 1) .NE.UJ GO TO 2 
RETURN 

2 NUPLO T = »E AL oE » 

Sx T HE T = • T R Ob . 

IE ( P L 0 1 ( 2 ) » b Q . 0 ) oKTHE T=. E AL SE • 

R b W IN 0 1 
NR EC= 0 
IP = 1 

WRITE to* 75) 

K E T UR N 
ENTRY OUK 

IF (NUPLO T) RETURN 

IE UCUTPLT) .ANU. (IP.tQ.LRbCMl) ) RE TURN 
IP = IP + 1 
RA( IE ) = k 

IF JELUTl 1) ,EQ. 0) GO TO 3 

V A ( IP ) = Su R T ( 3NGL l V2 ) ) 

R EA ( IP ) = jNGl ( RE ) 

3 IE ( PL 0T( 2 ) . t « » 0 ) GO TO 4 
XXA ( IP )=SNoL(Ul 1) ) 

Y Y A ( IP ) = SNbL ( Q( 2) ) 

ZZAi IP ) = 3NGL l Q ( 3) ) 

T HE TA A i IP )=ARCuS( LZAU P) / SNGL ( K ) ) *R T 0 

4 IE (PLOTi 3) .£(J. 0) GO Tu 8 
GO TO (o* 5) ,PC S 
RETURN TJ ORIGINAL VERSION 
CXI=- SINl SNuL J Q l 5) ) ) M SI N ( SNG L ( Q ( 4 ) ) ) 

SXI=SQKT( l.-uXI*CXI ) 

ETA=ATAN2(CUS( SNGL (Q l 5) ) ) , SI N ( S NGL < Q ( 5 ) ) ) *COS ( S NGL ( Q( 4 ) ) S ) 

GO TO 7 

6 CXI=COS( SNGL(Q( 5) ) ) 

SXt= 8 lN( SNGL ( Q( 5) ) ) 

E T A = SN GL ( w ( 4 ) ) 

7 CE TA = COS( E TA ) 

SETA=5 IN( t TA ) 

XAJ IP ) =SX I*oE TA 

Y A J IP ) =-iXl*CE TA 
Z A< IP ) =CX i 
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GAMMA A I IP ) = ARCUS l l SNGL I Oil) )*XA( IP) + S NGL (0(2) ) *YA< IP) 
1 + SNGL(d< 3) )*ZAiIP}| /(SNGL(K) ) ) * R TO 

8 IF (PLOT! 4).EU.O) GO TO 9 
IF (SKIhtT) GU TO 9 

THETA A I IP ) =ARCU 8 1 SNGL I U I 3 ) I /SNGL (Kl J *RTD 

9 IF (PLOT! 5) »E 0 . 0 ) GO TO 10 

PHIAI IP)=ATAN2( SNGL! U (2) 1 .SNGLIQU) I )*RTU 
IF I P H IAI IP ) • L T. 0. ) PHIA ( IP) =PHi All PH-360. 

10 IF UP.CT.LREC) RETURN 
IF INREC.GT.O) GO TO 12 
V I =VA 1 2) 

R DAB S (RA I 2 ) ) 

XXB =X X A( 2 ) 

YYB=Y YAI 2 ) 

ZZB =Z ZAI 2 ) 

XB=XA I 2) 

YB =YA I 2) 

Z B = Z A ( 2) 

00 11 11= 1 ,NP A 

11 0 A l 1, I 1 ) =UA I 2 * I 1) 

12 NREC=NREC +1 
WRITE (1) UA 
DO 13 Il=l,NPA 

13 UAI l, 1 IDUAILREC, ID 
IP = 1 

RETURN 
ENTRY OOP 

IF (NOPLDT) RETURN 
IF (NKEC.GT.O) GO TO 15 
R DAB SIRA ( 2) ) 

V I=VA I 2) 

XXB =X X A( 2 ) 

YYb=Y YAI 2 ) 

ZZB=ZZA( 2 ) 

XB= XA I 2) 

YB=YA I 2) 

ZB=ZA I 2) 

OU 14 11=1, NPA 

14 UAI 1, I I )=GA( 2 , I 1) 

15 IPL = I P - 1 

IF I IPl.GT. 1) GU TO 16 
IF I IPL.Eo.i) IPl =LREC 
IF I IPL. EG. 0) IPL=LRtC-l 
GU TU 17 

16 NREC =N REC + 1 

17 M T 1 = 1 

IF INREC.Eu.l) uO TO is 
M T 1 =2 

WRITE II) UA 

18 XXE =X X AI IPL) 

YYE=Y YAI IPL) 

ZZE=ZZA( IPL) 

X E =XA I IPL ) 

YE = YA I IPL ) 

ZE = ZA I IPL ) 

IF IPLUTI U.EU.O) GU TU 31 
CALL H EAO 

CALL L R SI ZE (0.0,10. ,5.0,10.) 

CALL LKMRGN I 1.0,0.5,0.5,0.51 
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CALL LrLEGN ( 36Hi UN-MOLECULE it PARA TI ON IN ANGi TROMS » 36 ,0, 3 . 84, 5 . I 
123, NEOP) 

lALL LRLEGN I 35HVELUC I TY SCALED TU INITIAL VALUE OF , 35 , 1, 0 . LOO , 5 * 3 
1 3 » N EO P ) 

call lrcn vt < vi *i .to, 3 , i v .,4 a i ,4) 

CALL LRLEGN ( IV, 11, 1,0.100,8. la , Nt UP ) 

CALL lRLEGN ( 7n CM/SEC , 7, I ,0. 100 , 9 . OS , NEUP) 

CALL L R A N G b I R I , 0® 0, 0.0,0. 01 
CALL LkGkIU 12,1,RI/10»,6«) 

GO TO U9,2l),MTl 

19 DU 20 11=1, IPL 

20 VAI ID =VA ( II ) /VI 

CALL LRCURV I K A , VA , I P L , 2 , S Y Mb OL , Nt uP) 

GO TU <:4 

21 REWIND 1 
N P = L R E C 

DO 23 IRE C = 1 » NKEC 
RtAD II) UA 

IF I IREC. tw.NktC) NP = I P L 
DO 22 I Jl = l * IMP 

22 VAI 1 1 ) =VA I I 1 ) /VI 

23 CALL LRCJRV ( R A » VA , NP » 2 » SY M8UL, NEUP ) 

24 CALL LRS1ZE 10.0,10. ,0.0, 5.0) 

CALL LRLEGN l 36H l UN-MULEC ULE SEPARATION IN ANGSTROMS, 36, Of 3.84,0.1 
12 5, NEOP ) 

CALL LRLlGN ( 32HR0TATI ONAL ENERGY SCALED TO NT = , 3 2 , 1 , 0 . 10 0 , 0 . 6, NE 
10P ) 

CALL L RCN V T (KT ,3,KTC ,4,11 ,4) 

CALL LRLE GN I K Tl , 1 1 , 1 , 0. I 00 , 3 . 1 8 ,NE UP) 

CALL LRLEGN I 2HE V , 2 , 1 , 0. I 00 ,4 . L 2 ,,NE OP ) 

REbCAL =K T * 1. 6fc- 4 
Gu TU I 23, 2 7), Mil 
23 DO 2b 11=1, IPL 

26 REAI 1 1 )=KEA( l 11 /kESCAL 

call LRCURV (KA*REA,I PL, 2 , SYMBOL ,NEOP ) 

GO TU 30 

2 7 REWIND 1 

NP =LR tC 

DO 29 IKEC = 1,NREC 
READ I l) wA 

IF UREC.eQ.NREC) NP=IPL 
DO 28 11=1, NP 

28 KEAI I 1 )=REAl 1 1) /RESCAL 

29 CAlL LRCURV I K A ,RE A ,NP ,2 , S Y Mb OL , NEOP ) 

30 CALL L RLE GN I 1H , 1 , 0 , . 5 , . 5 , E UP) 

31 IF (PLOT! 2) .EO.O) GO TO 48 
CALL HEAD 

CALL L RMR GN I 1.0,0. 3,0*5, 1.0) 

CALL L RSI ZE (0.0,5.0,5.0,10.) 
call l range I -r I , k I ,-r I ,r I j 
call l rgr id 12,2,ri,ri) 

CALL LRLE GN I oHX AXIS, 6, 0,2. 25, 5. 125, NEOP) 

CALL L RLE GN (6HY A XI S , 6 , 1 ,0 . 1 00 , 7. 2 5 , NEOP ) 

GO TO I 32 , 33) ,MT1 

32 CALL LKCOKV I X XA , Y YA , I PL , 2 , S Y Mb OL, NE OP i 
GO TU 35 

3 3 REWIND 1 

NP =LR EC 

DU 34 IRE C = 1 , NREC 
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RtALI (1) UA 

I F (iRtC.Eu.NRLC) NP = I PL 
34 CALL LRCURV { XX A » Y YA , NP , l , S V MBOL , NE UP 1 
33 CALL LRLAdL ( d t G , 4 , 0 » XX8 , V YB , NE OP) 

CALL LRLAdL IEND»4,0»XXE,YY£,NEUP) 

CALL L RSI Zt ( 3. 0. 10. ,5.0,10. ) 

CALL LRLtGN ( 6H Y A XI S , 6 , 0 * 7 .2 5 , 3 . 12 3 , NEOP » 

CALL LRLcGN I 6H Z AXIS, 6, 1,3. 100, 7. 25 , NEUP) 

GO TO ( 3 6 » 3 n » M T 1 

36 CALL LRCURV 1 Y Y A , Z ZA , I PL , 2 , S YMBOL , Mt OP » 

GO TU 39 

37 REwlNJ 1 
NP =LR EC 

DO 36 IRt C = 1 » NRtC 
READ { 1) wA 

IF i IRtL.Ew.NREC) NP=IPL 

38 CALL LRCURV ( YYA , ZZA , NP , 2 , SYMBOL NEOP ) 

39 CALL LRLAdL ( BEG , 4 , 0 , Y Yd , ZZd , NE OP) 

CALL LRLAtJL ( END ,4 » 0 , Y YE , ZZE , iME OP) 

CALL LRSIZE l 0. G» 5. U,0. 0,3. 0) 

CALL LRLtGN ( 6HZ A XI S , 6 , 0 ,2 . 2 5 , 0. 1 2 5 , NEOP ) 

CALL LRLEGN < 6H X A XI S , 6 , l ,0. L 00 ,2. 25 , NEOP) 

GO TO (40, 41), MU 

40 CALL LRCURV ( ZZ A , XXA , I PL , Z , S YMBOL , NE 0 P ) 

GO TO 43 

41 REWIND l 
NP =LR EC 

DO 42 IREC = L» NR EC 
READ ( 1) uA 

IE ( IREC.EO.NREC) NP=IPL 

42 CALL LRCURV ( ZZA , X XA , NP , 2 , S YMbQL ,NE UP ) 

43 CALL LRLAtJL ( BE G , 4 , 0 , Z ZB , XXfl , NE OP) 

CALL LRLAdL ( END , 4 , 0 , ZZE , XXE , IME OP) 

CALL LRSIZE ( 5.0,10. ,0.0,5. 0) 

CALL LRLtGN (36H10N MOLECULE SEPARATION IN ANGSTROMS, 36,0, 6.34, 0. I 
125, NEOP) 

CALL LRLEGN ( 26H PULAR ANGLfc IN DEGREES ,26 , 1 , 5 . 100 , 1 . 50 , NEOP ) 
CALL LRANGt(RI ,0.0,0. 0,180. ) 

CALL LRGR ID < 2 , 2, R I / 5. ,4 5. ) 

GO TO ( 44 , 43 ) , M T 1 

44 CALL LRCURV ( R A , THE TA A , 1 P L , 2 , SYMBOL , NE CP ) 

GO TU 47 

45 REWIND 1 
NP =LR EC 

DU 46 IRE C = 1 , NR EC 
READ (1) UA 

IF ( IREC.EO.NREC) NP=IPL 

46 CALL LRCURV ( R A , THE TA A , NP , 2 , S YMdOL , NEOP ) 

47 CALL LRLtGN ( 1H ,1,0, 5. 5,. 3,E0P> 

48 IF (PLOT! 3).EU.O) GO TO 65 
CALL HEAD 

CALL LRMRGN ( 1 . 0, 0. 3 , 0. 5 , 1 . 0) 

CALL LRSIZE (0.0,5.0,5.0,10.) 

CALL L RANGE ( - 1 .0 , 1 . 0 , - 1 . 0 , l . 0) 

CALL LRGR1D ( 2 , 2, 1.0, 1 . 0) 

CALL LRLEGN < 6HX A XI S , 6 , 0 ,2 . 2 5 , 5 . 1 2 5 , NEOP ) 

CALL LRLtGN ( 6HY AXIS, 6, 1,0.1 00, 7. 25, NEOP) 

GO TO 149 , 50) ,M T1 
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49 CALL LKCURV l AA , YA , I P L» 2 » SY Mti UL , NE OP) 

GU IU 52 

30 K LW IN 0 1 
NP =LR LC 

Du 51 IRfcC = 1 » NK EC 
READ (II 4 A 

IF ( IREc.fcU.NRfcC) NP = 1 PL 

31 CALL LkCURV ( XA,YA,NP,2 » SYMBUL,NEOP) 

52 CALL LRlA bL (BEG,4,0,XB,YB,NEOP) 

CALL LRLAdL ( END , 4 » 0 » XE , YE » NE OP ) 

CALL LRS1ZE ( 5.0, 10. ,5.0,10.) 

CALL L RLE GN (oHY A XI 3 , 6 » 0 » 7 . 2 5 , 3. 1 2 5 » NEUP ) 

CALL LRLEGN (GHZ AXlS»6,l,5.l 00 , 7 . 25 , NEUP ) 

GO TO (33, 54), MTI 

53 CALL LKCURV ( YA » ZA , l P L» 2 , SYMb UL , NE UP ) 

GU TO 5o 

54 REWIND 1 
NP =LR EC 

DO 35 IREC =1 , NREC 
READ ( 1) UA 

IE ( IREC. EO. NREC ) NP=IPL 
35 CALL LKCURV ( YA ,ZA , NP , 2 , S YMB UL, NEOP) 

56 CALL LRLABL ( BEG , 4 , 0 , YB , ZB , NE OP ) 

CALL LRLABL < END , 4 , 0 , YE , ZE , NE OP ) 

CALL LRSIZE 10.0,5.0,0.0,5.0) 

CALL L RLE GN (6HZ AXI S , 6 , 0 ,2 . 2 5 , 0 . I 25 , NEUP ) 

CALL LRlEGN ( 6HX A XI 6 , 6 , 1 , 0. 1 00 , 2. 25 , NEOP ) 

GU TU (57, 36), MTI 

57 CALL LRCJRV < Za ,XA , IP L,2 , SYMBOL, NfcUP) 

GO TO 60 

5 8 R tw IN J 1 

NP=LRtC 

DU 39 1R E C = 1 » N R Eu 
READ ( 1) wA 

IE ( IREC.Eu.NRLC) NP = IPL 
39 CALL LKCURV l ZA , XA , NP , 2 , S YMB UL , NEUP ) 

60 CALL LRLABL ( bL G , 4 , 0 , Zb , XB , Nt UP ) 

CALL LRLABl ( END , 4 , 0, ZE , XE , NEOP ) 

CALL LRS1ZE ( 3. 0 , 1 0. , 0 . 0, 5 . 0) 

CALL L RLE uN ( 3ohI UN-MULEC ULE SEPARATION IN ANGSTROMS, 36,0, 6.34,0.1 
123, NEOP) 

CALL LRLcGN < 2BHuR 1 EN I A Tl UN ANGLE IN DEGREES , 2 B , 1 , 5 . 100 , l. 50, NEOP ) 
CALL LRANGEU l , 0.0,0. 0, L BO. ) 

CALL LRGkIO ( 2, 2,RI/3. ,45. ) 

GO TO { 6 L , o2 ) , M Tl 

61 CALL LRCJRV ( R A , GAMMA A , l P L , 2 , SYMBOL , NEOP) 

GU TO 64 

6 2 REWIND 1 

NP =LR EC 

DU 63 IRE C = l » NREC 
R EAO ( l i 4A 

IE ( IREC.EU.NRLC) N P = I PL 
63 CALL LKCURV ( R A , uAMMA A ,NP , 2 , S YMBCL , NE OP) 

6 H CALL LRLEGN ( 1H , l , 0, 3 . 5 , . 5 , E OP ) 
o 5 IE (PLUK 4J.CU.0) GU TU 70 
CALL HcAU 

CALL CRSIZE ( D. ,10. ,0. ,10. ) 

CAlL LkMkGN (1. 0,0.3, 1.0, 0.3) 

CALL LRCMSZ (2) 
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CALL LrtLEuN ( 36HiuN-MULECULE SEPARATI ON IN ANGST ROMS, 3t> , 0» 3 . 14, 0 .2 
100* NEUP ) 

CALi_ L RLE GiN { 22HP0LAK ANGLE IN UEGkc E S ,22 , 1 » 0 . 12 , 3 » 80 » NEGP ) 

CALL LRANGEIR i ,0.0,0. 0,180. ) 

CALL L K Gk 1l) {2,2?RI/t>«,45.) 

GO Tu I 60 , 67 ) » M T 1 

66 CALL LRCUKV ( k A , THE TAA , l PL , 2 , SYMBOL , NE CP) 

GU TO 69 

67 REWIND 1 
NP =LR EC 

00 68 IKE C = 1 , NkEC 
READ (1) DA 

IF i iREC. tw.NREC) NP=IPL 

68 CALL LKCUkV (RA»ThETAA,NP,2 , SYMBOL, NEOP) 

69 CALL LRLtGN l 1H ,1,0, *5, .2, EOP) 

70 IP I PL 0 U 5) .Eu.O) GU ID 100 

CAlL head 

call LRSl Zt ( 0. ,10. , 0 . , LG. ) 

CALL LRMRGN 11.0,0.5,1.0,0.5) 

CALL LRCH 8 Z (2) 

CALL LRLtGN 1 36HI UiN-MULEC ULt SEPARATION I N ANGST ROMS, 36,0, 3. 14, 0 .2 
100, NEUP ) 

CALL LRLtGN ( 26HA ZI MU THA L ANGLE IN DEG REES, 26 ,1,0. 12,3 .AO, NEQP) 
CALL LRANGtlRI ,0.0,0. U,3 oO.) 

CALL LRoklU I 2 , 2 *k I /5 • ,45 . ) 

GO TO (71,72)»MTL 

71 CALL LRCUR V ( R A , PH I A , I P L , 2 , S YM 8 UL , N EOP ) 

GO TO <4 

72 R Ew IN D l 
NP =LR EC 

DO 73 1 RL C = i ,NkEC 
READ ( 1) UA 

IF i IREC.EQ.NREC) NP = IPL 

73 CALL LRCURV i R A ,PH I A , NP , 2 , S YMBOL ,NEOP ) 

74 CALL L RLE GN I 1H ,1,0, .3*. 2, EOP) 

100 CONTINUE 

RETURN 

75 FORMAT (24H+IHIS CASE TO BE PLOTTED) 

END 


$I8PTC HEAD. 

SUBROUTINE HEAD 

COMMON /ODC/ RE,V2,SXI ,R, T , KT , I D 1 , I D2 , U < 12 ) , CLRP , RU N( 3 ) , PL QT ( 5 ) 
1 , C U I P L T 

DOUBLE PRECISION KE,V2 ,SXI ,R , T,KT, u 
REAL CASE I 2) ,NEUP ,TYPE 18) 

DA f A N EOP / 0. / , 

1 TYPE/ A 8HCAP TORE kEPULSION OUT OF TIME TOO SMALL / 

INTEGER CLRP, PLOT 
CALL L RSI ZtIU.0,10. ,0. 0,10. ) 

CALL L RCN V T I ID 1 , 1 ,CagE 1 1 ) ,1,6,0) 

CALL LRCNVTl 1U2,1,CASE(2) ,1,6,0) 

CALL LKCHSZI2) 
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INDEX = 2*CCLRP-1) * 1 

CALL LRLEGNI 3HR UN » 3* 0* 1. 50, 9. 875 ,NEOP) 

CALL LRLEGNI RUN, 1 8 , 0,2 . 00 , 9. 875 ,NEOPI 
CALL LRLEGNI 4HC ASE , 4, 0 ,4. 80, 9.875 , NEOPJ 
CALL LRLEGNICASE, 12,0,5. 425 ,9. 875, NE OP) 

CALL LRLEGNI TYPE! INDEX) » 12 , 0 » 7. 125 , 9. 875 , NEOP ) 
CALL LRCHSZI 1) 

RETURN 

END 


SIBFTC PLTCM. 


SPECIAL MUVIE DECK, SIMPLIFIED 
SUBROUTINE DOS 
LOGICAL UN 
LOGICAL CUTPLT 
LOGICAL UIPOLE 

REAL XAU500J. YAI1500) ,XXA(1500» , Y YA ( 1500 ) , T AI 1500) ,NEOP, 

1 DUMXI 2) , D UMY ( 2), ZR 1 1500) 

DOUBLE PRECISION E ROC, RM , V, ETRC ,E PLC , MUE , XC A.XCB ,XC , CD1 
DOUBLE PRECISION 0, V2,R,RE ,T,KT,SXX 
INTEGER PCS, CLRP, PLOT 

COMMON FOR POTENTIAL SPHERE RADIUS 
COMMON /CIRCLE/ RADI US,RV IEW ,0N 

COMMON FOR COORD SYSTEM SWI TCH 
COMMON/PC SC/ PCS 

COMMON FOR ALL MAJOR VARIABLES 

COMMON /DOC/ RE ,V2.SXX,R, T,KT, I 01,102,01 12) , CLR P, RU N< 3) , PLOT ( 5 ) 
1 .CUTPLT 

COMMON /J AC KEN/ ERliC , RM , V.ETRC , EPLC , MUE , XCA , XCB»XC, JCD, CD1 
D IMENSIUN XAXI SXI2) , YA X I S X ( 2 ) , YAXl SYI2) ,XAXISY<2) 

DATA XAXI SX/-10..+ 10./, YAXISX/-10. ,-10./ 

DATA XAXI SY/-10.,-10./, YAXI SY/-10. ,H0. / 

DATA AX,BX,CX,DX,EX /3*2. 0,5. 0,8. 0/ 

DATA AY,BY,CY,DY,EY /2.0, 5. 0, 8. 0,2. 0,2 . 0/ 

DATA LKECM1/1500/, EUP/1./, NEOP/0. / 

DUMMY ENTRY TO START PLOTTING 
MAKES MOVIES IN X — Y PLANE NOW 

IP = 0 
ON = * TRUE. 

DIPOLE = .TRUE. 

SAVE COLLISION RADIUS 
ft ADMAX=RAD I US 

I PI MUE .fcj .0 .DO ) DIPOLE = .FALSE. 

WR ITEI 6,101) 

101 FORMAT k 24H+M0VIE MADE IN CM SYSTEM 1 
R ETURN 

ENTRY TO RECORD VARIABLES 
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ENTRY DOR 

I FI IP «GE®LRECM 1 ) RETURN 
IP » IP * 1 

STORE SCALED TIME 
CLOCK PERIOD IS 10##4-i3I SEC 
T Al I P I = AMODI SNGL IT) *I0« 3 **62 6319 

STORE ION POSITION 

CENTER OF MASS COORDS 
X = 04 II 
Y = Qi?) 

PHI = ATAN2 4 Y* X J 
XXAIIPI = 0*5# SNGL I R 3 #CGS I PHI I 
YYAUPI = 0®5*SNGLlR)*SINlPHI) 

Z = 0* 5*0 I 31 

IF4ABS(Z)*GT«io«) GO TO 2 
SPECIAL RADIUS VARIATION 
RADIUS = RADMAX*4G® 5*Z/40®I 
GO TO 3 

2 RADIUS = 0*0 

3 ZR< IP) = RADIUS 

STORE DIPOLE POSITION 
I FI *NOT*D I POLE ) RETURN 
GO TO 4 6* 5) * PCS 

5 CXI -SINISNGLIQI 5))J*SI N4 SN3 LI 044) II 
SXI = SORT* 1®0-CXI**2) 

ETA = ATA N2 4 COS I SNGLIO 4 51 I I * SI N4S NGL I OC 5 ) I 1 *COS 4 SNGL 1044) I ) I 
GO TO 7 

6 CXI = COS C SNGL 4 01 51 I ) 

SXI = SIN 4 SNGL 4 015) ) ) 

EIA = 0441 

7 CETA * CD SI E TA I 
SETA = SIN4ETA) 

NOTE SCALING OF DIPOLE LENGTH 
*57735 IS I® 0/ SORT! 3* I 
SCAL ER= *57 73 5*R ADI US 

NO It PHASE SHIFT BECAUSE ETA AHEAD BY PI/2 
XAIIPI = SX I #SE T A# SCALER - XXA4IPI 
YAIIP1 = ~SXI*CETA* SCALER - YYA 4IPJ 
RETURN 

ENTRY TO MAKE MOVIE 

ENTRY OOP 
£ CLOSE IN VIEW 

R VIEW = 10* 

CALL LRMON 
CALL LREQN 

CALL LK SIZE I 0®Q*I0® *0®Q*10® ) 

CALL LRMRGN4 <,75*® 75*®75*» 75) 

CALL BLANK 
CALL TEST 
CALL BLANK 
CALL TEST 
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CALL L KMRGN { 2. 350,2. 0,2. 233,2. 0) 

CALL HEADM 
CALL LRCHSZ ! 2) 

CALL LRANGEl— KVIEW*RV1EW*— RVIEw.RVIEW) 

CALL LKGRIDI 1,1,0., 0.) 

BEGIN MAIN MOVIE LOOP 

DO 80 1=1, IP 
I SAVE = I 

MAKE SURE WE HIT OUTPUT 
CALL T IMLFTS TLEFT) 

I f < TLEET.LT. 900. ) RETURN 

LABEL THE PERIMETER (IN ANGSTROMS) 

CALL LKLEGN ( 3H-10. 3 . 0. AX . A Y , NEOP ) 

CALL L RLE GN I 1H0 » 1 , 0,8X ,B Y,NEUP ) 

CALL L RLE GN( 1H0 ,1 , 0,DX,0Y,NEUP) 

CALL L RLEGN ( 2H 1 0 *2*0, CX,CY» NEOP) 

CALL L RLEGN 4 2H 10 .2, 0,EX,EY,N£ OP) 

DRAW THE AXES 

CALL LRCURVt XA XI SX» YAX I SX *2*2 , SYMBOL » NEOP) 
CALL LRCUR V ( XAXISY,YAXISY ,2*2* SYMBOL* NEOP) 


PLOT ION ONLY IF NEAR ENOUGH 

X TEST = XXA I I) 

Y TEST = YYA (I) 

I FI ABS4 XTEST 1.GE.RVIEW .OR. ABS(YTEST) .GE.RVIEW) GO TO 78 
CALL IRC UK VI XTE ST* YTEST* 1 *3 ,1H + *NE0P) 

GET RADIUS FOR DOTS 
RADIUS = 7RI IS 
IFIRAOIUS.EQ.O. ) GO TO 7B 
CALL DO ISl XTEST, YTEST) 

78 CONTINUE 

PLOT THE CLUCK 

X CLOCK = 0. 5*S INI TA4 I ) ) * 8. 

Y CLOCK = 0 . 5#C US4 T All) ) + 8. 

CALL LRLEGNI 1H*. 1, 0,8. ,8. ,NEOP ) 

CALL LKLEGN ( 7H$GNG$GE ,7,0, XCLOCK*YCLOCK, NEOP) 


PLOT THE DIPOLE 
I F( DIPOLE) GO TO 20 

NO DIPOLE — MUE=0 
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€ PLOT CENTER 

CALL LRCUR V(— XTE ST«— YTEST *1 ?3 * 1HG»N£QP) 

CALL OOTSI-XTESTe-YTEST) 

GO IQ 75 
C 
£ 

£ NEW COMPLICATIONS 

£ XTEST* Y IE ST ARE COORDS OF NEG END 

C X $ Y ARE COORDS OF POS END 

C 

20 CONTINUE 

X IE ST = XA (IJ 

Y TEST - YA I I ) 

X - -XTEST - 2« G^XXA I 1 1 

Y = -YTEST - 2« Q#Y YA III 
I BOTH = 0 

I FI ABSI XTE ST UGE^RVIEW » OR* ABSI YTEST) . GE. RVI EW) GO TO 31 
I BOTH = I BOTH * 1 
£ NEG END ON 

CALL LRCURVi XTE ST* YTEST* 1*3 *1H-*NEUP) 

31 IFI ABSI X ) *G£ ®R VIEW *GR* ABSIYI.GE.KVIEW) GO TO 32 
I BOTH = 1 BOTH * 1 

C POS END ON 

CALL LRCURV4X*Y*1*3*1H**NE0PI 
CALL LRLABL ( 7H$GNW$GF t 7 #0 , X* Y» NEUP * 

32 IFI IBiITH.LE.lJ GO TO ?5 
IFI RADI US «E 0*0* ) GO TO 75 
RADIUS = RAGHAX - RADIUS 

I FI X TE ST® L T*> X) GO TO 76 
DUMXI1I = X 
0 UM X ( 2 ) = XTEST 
DUMYIU = Y 
DUMYI2I = Y TEST 
GO TO 77 

76 DUMXI II » XTEST 
DUMXI 2 I = X 
DUMYI II = YTEST 
DUMYI2J = Y 

77 CONTINUE 

CALL DOTS I-XXAII)®— YYAIIt ) 

CALL LRCURVI DUM X*DUMY *2* 2 * SYM3 GL* NEOP ) 

7 5 CONTINUE 

CALL LRCURVI 0* 0*0* 0* U 3* 1H *£QP) 

BO CONTINUE 
CALL BLANK 
CALL LREOFF 
CALL LRMOFF 
CALL L RCHSZ III 
RETURN 
END 


45 



SIBFTC HE: A DM 


SUBROUTINE HFADM 
C HEADER F OH MOVIES 

I NT EGER CLRP .PLOT 
REAL CASE 12) ,TYPE( 8) 

OUUBLfi PRECISION RE. VV.SXl ,R,T,KT,Q 

COMMON /DOC / RE . VV. SX1 »R»T»KT,ID1, 102,01 12) , CLRP, RUNl 3 ) , 

* PLOT I 5 ) tCUTPL T 

DATA TYPE/48HCAP CURE REPULSION OUT OF TIME TOO SMALL / 

CALL LKGHSZI4) 

CALL l.RCN V T ( ID1 ♦ l.CASE ( 1) ,1,6.0) 

CALL LKCNVT 1 102. 1 .CASE (2) *1 ,6,0) 

INDEX = 2*1 CLRP— 1) * i 
DO 10 1=1.20 

CALL L RLE GN < 3HR UN«3»U*1. 0,9.0, 0.0) 

CALL L RLE ON 1 RUN. 18, 0,2. 5 ,9. 0,0.0) 

CALL L RLE GN ( 4HC ASE , 4, 0 ,1. 0,7. 0,0.0) 

C ALL LKLEGM CASE, 12, 0*2. 6,7.0, 0.0) 

CALL LRLEGN1 TYPE { INDEX) ,12,0,1.0,6.0,1.0) 

10 CUNT INUE 
RETURN! 

END 


SIBFTC DOTS. 

SUBROUTINE DOTS ( XCENTR , YCENTR) 

LOGICAL ON 

REAL P T SI 6) » PT 14) 

COMMON/CIRCLE/ K,RM,ON 

DATA PT/i.O, .92388, .70711, .38268/, PTS(5)/0.0/ 
C 

X = XCENTR 
Y = YCENTR 

I F( .NOT .UN ) GO TO 20 

C TEST CORNERS 

C = X+R 

IF1C.GE.RMJ GO TO 50 
C = X-R 

IFIC.LE.-RMI GO TO 50 
C = Y*R 

I F1C.GE.RM) GO TO 60 
C = Y-R 

IE1C.LE.-RM) GO TO 50 
C CORNERS ON THE PLOT 

20 CONTINUE 

C THIS IS FOR VARIABLE RADIUS 

DO 25 1=1,4 

25 PIS( I) = PT< I)*R 
30 DO 40 1=1,5 
J = 6-1 

CALL LRCOKVl +PTS ( J ) + X , *P T SO ) + Y , 1 ,3 *1 H* , NEOP ) 
CALL LKCUR V ( +P TS ( J ) +X ,-PT S (I) * Y , 1 , 3 , l H* , N60P ) 
CALL LRCURVl-P TS1 J)*X,-PTS1 I )*Y,1 ,3 ,1H*, NEOP) 
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CALL LRCURVI-PTSI JH-X.+PTSI I ) + Y, 1 , 3 »1H* » NEOP ) 
40 CONTINUE 
50 RETURN 
END 


$ IBFTC BLANK. 

SUBROUTINE BLANK 

REAL XU2 )«X2< 21 , X3I2) »YU2),Y2t2) ,Y312) ,Y4 (2) 

DATA XI /-l.C.l.O/, X2 /l. 0.1.0/, X3 /-1.0.-1.0/, Yl/-. 766, -.766/ 
DATA Y 2 /- . 766« .766/, Y3 /.766,. 766/, Y4 /. 766, -.766/ 

MAKE BLANK FILM 
00 10 1=1,80 

CALL LRCURV< 0.0,0.0,0.1,1HX,1.0) 

10 CONTINUE 
R ETURN 

ENTRY TEST 

MAKE TEST PATTERN 
CALL LRGRiDI 1, 1,0. .0. I 
CALL L RANGE (-1.0, 1.0. -1.0, 1.0) 

CALL LRCURV(X1, Y 1,2.2. SYM, 0.0) 

CALL LRCURV (XI, Y3 .2, 2, SYM, 0.0) 

CALL LRCURV! X1.Y2.2.2 , SYM. 0.0) 

CALL L RCUR V ( XI ,Y4, 2,2. SYM, 0.0) 

CALL LRCURV1 X3. Y2.2.2.SYM.0.0) 

CALL LRCURV ( X2*Y2<2,2 .SYM.1.0) 

RETURN 
END 


IBFTC PLTOUM 


DUMMY ENTRIES WITH PURPOSES AND ARCS EXPLAINED 
MORE INFO IN NASA TM X-1866I 1969) 


***** MAIN ROUTINE. DRAWS A CURVE T 

SUBROUTINE LRCURV (X.Y .NPTS.NTYPE .SYMBOL, EOP) 

X— ARRAY OF ABCISSA VALUES 

Y — ARRAY OF ORDINATE VALUES 

NP T S — NUMBER OF THESE VALUES 

NT YPE — TYPE OF PLOT. .. 1 =POI NTS ,2 = LI NES » 

3=SYMBQLS .4= SPECIAL SYMBOLS 
SYMBOL-- THE PLOTTING SYMBOL 

EOP — END OF PLOT SMI TCH « » . 0 . 0= NOT END, 1.0= END, ADVANCE 



c 

t 

c 

c 

c 

c 

C 

L 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

(,***** 

C 

c 

l 

C 

C 

c 

c 

c 

c 

c 

c 

c 

C 

c 

c 

c 

c#£* 

c 

c 

c 

£***£* 

c 

c 

c 

c 

u 

c 

c 


PR Ii\T^> UJT L t G F ND S 


EN T RY lkL t GN ( C HAR S ? N # NOR 1 E N f X ? Yj EUP ) 


CHARS" 

N — 

NU R It H — 
X — 

Y — 

EJP — 


ARRAY l JF CHARAC FEKsIBCD) TU BE PR I NT ED (6 PER tiUKD) 
NUmGEk JF CHARACTERS* ®®£®G® i2=CHARS DIM® BY 2 
JkIENTA JIUN. * * 0=HURI ZON T AL , i = VERT i CAL PRINT IN G 
X CJORU OF FIRSI CHARACTER ? IN Abb® UNITS 0® — 10® 

Y CU UK J uF FIRST CHARACTER, IN ABS® UNITS 0.— 10. 
AS AbUVE 


DEFINtS USER 8 S UNITS FuR PLOTTING 
ENTRY LkANGL ( XMIN ? XMA X ? YM I N * Y MA X ) 

XM IN-- rtllM X VALUE? USER'S UNITS 

XMAX — MAX X VALUE? USER 8 S UNITS 

YM IN-- MIN Y VALUE? USER'S UNITS 

Y M AX — MAX Y Value? USER'S UNITS 


LIRE LRLLUN? BUT WORKS IN USER'S UNI TS 
ENTRY L R L A B L ( C H A R S » N ? H OK i t N ? X ? Y ? E 0 P ) 

ARCS SAME AS FOR uKLEGN? EXCEPT X AND Y NOW IN JSER 9 S UNITS 


DEFINES MARGINS FOR This PLOT 
ENTRY LkMRGNI XLeF T * X* 1 GH T * YBO TM ? YT OP ) 

XL L FT — LEFT MARGIN IN AGS® PLOTTING UNITS 

XKIGHI— RIGHT MARGIN IN AbS® PLOTTING UNITS 

YBO TM-~ BOTTOM MAkGIN IN ABS® PLOTTING UNITS 

Y T U P — TUP MARGIN IN ABS* PLOTTING JIM ITS 


defines grid fur this plot 
ENTRY LRokIlH IX?IY,DX?UY) 

IX — CODE FUR VERTICAL GRID LINE 0= S T ANOAR D? 

1 =0 X GIVES NO® OF LINES 

2 =D X GIVES INTERVALS BETWEEN LINES 

3 =D X GIVES NO® OF TICK MARKS 

4=UX GIVES INTERVAL BETWEEN TICK MARKS 
I Y — Cuuc FOR HORIZONTAL GRID LINES? USED AS IX IS USED 

DX — NUMBER OR INTERVAL? VERTICAL LINES 

DY — NUMBER OR INTERVAL? HORIZONTAL LINES 


DEFINES PIECE UF t*k A ME FOR THIS PLOT 
ENTRY LRSIZfcl XL EFT ?XRiGHf ?YBOTM?YTUP) 

ar^s as for lrmrgn 


CONVERTS NUMBER TU BCD FOR PRINTING 
ENTRY LRCwVn VALUE , NT YPe ,CH ARS? NF ORM? N ?M) 

value — NUMBER Tu BE CONVERTED* TYPE IS NT YPE 

NT YPE — ® ® ® L — VAL UE IS l fslTEGER ( F I XEu) * 3= V ALUC IS FLOAT ING 

CHAR S — ARRAY Tu HOLD CHARACTERS INTO WHICH VALUE IS CGNV * 

NFJRM — FORMAT C UDE . . . 1 = I ? 3 =F ? 4= L 

N — TuTAL NU® UF CHARACTERS 

M — CHARACTERS TO RIGHT OF DECIMAL* = 0 FUR NFU RM= i 
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o o r r o o or r o r c o o 


***** TURN* tXTtNDE J MUut l UiE G rtHUlt FRAME ) ON AND UFF 
EN TRY LRE ON 
ENTRY LrEGFF 

***** turns movie mjde i ruta te g vo degrees) on and off 
entry lkmon 
EN T t\Y lRMUFF 


***** GeTG CHARACTcR SIZE 
ENTrY LRu HSZ INSIZE) 

NSIZE — DEFINES CHARACTER GIZt, 0 = L L T PLUfTING ROUT. PICK, 

I = T I N Y , Z = G MAL L » 3= MED* 4=BIG 


ENTRY TIMLFTI TIME ) 

***** reads INTERNAL SYSTEM C LUC K , RE TURNS TIME LEFT TU END OF RUN 
TIME— TIME REMAINING, IN PuLiES. . . ONE PJL,SE=i/oO SECOND 


ENTRY AREKR(TEXT) 

C**#** FOR TKAN-CALLtD ARITHMETIC ERROR ROUTINE IN OUR SYSTEM, 
C****# IRE ATE U LIKE LIBRARY ARlTrt ERRORS. ..CALLED ONLY BY UAKCOS 

C TEXT— MESSAGE TJ BE PRINTED WHENEVER CALL IS CXECUTED 

C TLXI MUST END IN $ 

C 

c 

K E I UR N 
END 


CONCLUDING REMARKS 

A FORTRAN IV program has been developed to numerically study ion - polar- 
molecule collisions. Results are obtained in the form of time history plots and motion 
pictures, as well as collision statistics. This report describes, for potential users, the 
structure and use of the program. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, July 29, 1970, 

129-02. 
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APPENDIX - PROGRAM FOR GENERATING INITIAL CONDITIONS 


The main program generates batches of random input conditions in the format ex- 
pected by the collisions program (on the Q cards). A MAP listing of the random number 
generator RAND/ SAND is also included. 


Input 


Two kinds of NAMELIST are used, one in the main program, the other in SUB- 
ROUTINE GENARY. They contain 


Block 

/INI/ 


Variable Variable in 
main program 


Meaning 


NC 

PSW 

BETA1 

BETA2 


number of cases of random coordinates for this 
set of values 

= 0, do not punch cards 

= 1, punch cards 

= 1/KT^; where KT ^ is the most probable 
energy, first degree of freedom 

= I/KT2; where KTg is the most probable 
energy, second degree of freedom 

BETA2 = 0 for linear molecules 


/GENER/ ATE MIA 

MIB 

RA 

VA 

BA 


first moment of inertia, 1^ 
second moment of inertia, Ig 
initial separation, r^ 
initial velocity 
impact parameter, b 


Subroutine GENARY 

This subroutine reads the array ATE (42) and scans it looking for the form limit 
increment j, limitg, incrementg, . . . and then generates values of uniform increment; 
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between limit^ and limit i+ j. It can be used to vary any or all of the parameters MIA, 
MIB, RA, YA, BA, in order to study the effect of these parameters on the trajectories. 


Main Program 

This program generates, in a series of nested loops, the coordinates and time 
derivatives, in sets of NC values for each set of parameter values. 


Output 

The main program prints the parameter values (in NAMELIST/ OUT 1/) and the 
random coordinates. Optionally, the coordinates are punched onto cards, suitable for 
input to the collision program, and the number of cards is printed (in NAMELIST/OUT2/). 
The coordinates are stored as follows: 


Mathematical name 

Stored in generator 
program in- 

Stored in collision 
program in- 

R 

(not used) 

Q(l) 

e 

A(n, 1) 

Q(2) 

<p 

A(n, 2) 

Q(3) 

v 

A(n, 3) 

Q(4) 


A(n, 4) 

Q(5) 

* 

A(n, 9) 

Q(ll) 

R 

DR 

Q(6) 

e 

AA(n, 1) 

Q(7) 

<p 

AA(n, 2) 

Q(8) 

V 

AA(n, 3) 

Q(9) 

k 

AA(n, 4) 

Q(10) 


AA(n, 5) 

Q(12) 
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Listings 


The following is a complete listing of the main program, subroutine GENARY, and 
subroutine RAND/ SAND: 


3.16PTC RANIT 

REAL A (999, 10 ) , RA( 1U J ,MI A { 10) »MIB( 10) , VA< 20) , BA( 133) , AA(999, 13 ), 
1 4 1 1 , M 1 2 

INTEGER PS'ni 
DATA TwUP 1/6.28 31853/ 
wKl TE< 6, 1001) 

loot Format ( ih i ) 

GALL SANG ( STORE) 

NAM tL IST/INi/ NC,PSW,8ETAl,BETA2 
5 REALM 3, IN 1) 
rtR l TE l 6, i Ni) 

P Srt-2-PSW 
GJ 10 1C=I,NC 
CALL RANG(RPHI ) 

R THETm = . 5 
CALL RANG ( RE T A ) 

CALL RANG(KXI) 

CAlL RANG ( RPS I 1 ) 

C A<_ L KANG ( RPS 1 2 ) 

PSi i=T*uP I*RPSI 1 
P SI 2= TwQP I*RP SI2 
L THETA 

A( 1C, 1 )=ARCCGS( 1.-2. SR THETA) 

C PHI 

A ( IC, 2 ) = TNGPI *RPHI 
C ETA 

A ( IC, 3 ) = r*OPl*R£rA 
C XI 

A ( IC, 4)=AKCCjS( 1.-2.SRXI) 

A( IC, 5 ) = S INI PSl i) 

A ( lC,6)=CUS(PSi 1) / ( 2. * SQR T ( R THE TAM 1 . -RTHET A) } ) 

A ( IC, n=CUS(PSI2)/{2.*SURT<RXI*(l.-RXim 
A ( IC, 9 )=0 . 

I F ( 8ETA2.EU.0. ) GD TG 10 
CALL RANG (KPS I 3) 

C PS I 

A ( IC, 9 ) = TWGPI*RPSI3 
10 A ( IC, d )=S INtPSI 2) 

CALL GEN4K Y( M I A ,NM l ,10,1 7 H MOMENT UF INERTIA,17) 

CALl GENAkYI MIB >\fM2» 10,24H$ECOND MUMENT OF INERTIA ,24) 

CALL GtNARY(RA,NR,lO,9HINI TIAL R,9) 

CALL u>tNARY( VA,NV*20,16H INITIAL VELOCITY, 16) 

CAl. L GENARY( 3A,Ntt» 100, 16. -I IMPACT PArAME TER , 16 ) 

NCARG S=NC *NMI *NR*NV*Nb*i\IM2 
NAM EL I iT/ OUT 2/NCAR03 
1 F ( PS W. LG. DrtR I TE (6, JUT2) 

I G= 0 

GU 30 1M1 =1, NM I 
Ml 1=MI A( l MI ) 
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DU 30 IM2=1,NM2 
M I 2 =M I » ( I M 2 I 
00 20 1C=1,NC 
CALL KAN) ( RR E ) 

Kb l=-ALJo ( KRE ) /lit: TA 1 
F 2=SwR l'< 2.* Kb i/MIAC IMI)) 

F 3-0. 

AA ( IC 9 5 ) = 0 . 

C bTA DOT 

A A l IC , 3)=F2*Al IC, n 
C XI OUI 

AA I 1C, A) = F2*A( IC,B) 

i F( 8feTA2.cg.o . ) jU ru 20 
C Pi I 03 T 

A A ( IC , 5 ) = SJRK AE1/M12) - AA ( I C , 3 ) * ( l . -2 . *RX I) 

20 CONTINUE 

00 50 IK= 1 » NR 
K=KA( IK ) 

R 2 =R*K 

00 80 IV=1,NV 
V =VA( IV) 

UU 30 IB = 1 , NB 
B = UA( IB) 

P 1= VA ( I V ) $BA I IB)/R2 
C K OUT 

Ok = -SOKT< l.-BAI Io)*t3A(Iri)/K2)*VA<lV> 

00 JO 1 C = i , N C 

C THETA DOT 

AA l IC, i) = Rn( I C , 5 ) 

c phi our 

30 Ail IC , 2 ) - F 1*4 I IC , 6 ) 

C 

NAM EL I Si/UUTl/M 1 1,MI2 ,KEL,RE2 , V , B , F l , F2 ,F J 
wRI TE ( 6 , 3 UTi ) 

1 0— 10 + 1 

kRITEI 6,1002) < ( A I IC, 1 ) ,1 = 1 ,4) , 0 R , ( A A ( I C , I ) , l = 1 , 4 ) , AA( IC, 5 ) « A( IC, 9 
1 ) , 1 G, I C , IC = 1 , NC ) 

100 2 FORMAT { 1H0, IX, 3HTHET4 , 4X , 3H PH I , 4X, 3HE T A, 5 X , 2 HX I , 10 X ,2H0R, 3 X , 
L6H0THETA, 8 X , 4H0PH I , 8 X , 4HD E TA , 9X » 3HDXI ,3X»4HDPSI ,4X,3HPSI/ 

2( OP 4F l .3, lPt>E12.3»0PF7.3*4X»2l5) ) 

UU TO (AO, 80), PS* 

c 

40 PUNCH 1003, ( (A( IC , I ) ,1=1,4) , OK, (AA( IC,I) ,1=1,41 , AA(IC»5), A(IC,9) 

1 , I U , IC, I- =i* NC ) 

1003 FORMAT (OP 4F 8 . 3, 1? 6E1 0 . 3/OPF 5. 3»4X» 21 3) 

PUNCH 1004 

1004 FORMAT! lrtl/ ) 

50 CONTINUE 

00 TO 3 
END 


53 



$ I8FTC GENAR 


SUBROUTINE GENARY { A, L S N» LABEL, NL) 

REAL AIN) , ATE { 42) »L43 ELU ) 

DO 5 1=1,42 
5 ATEI I )=0. 

NAMELIST /GENER/ ATE 
NW=NL / 6 

IF! MOD(NL,6).NE.O) NW*NW+1 
WR I TE 1 6, 1 002 ) ( LABEL! I ) ,1=1 »NW) 

1002 FORMAT UHL, 47X.6A6) 

READ! 5, GENER I 

DO 10 M=1 ? 21 

IF(AT£(2*M).EQ.0») GO TO 20 
10 CONTINUE 

CALL ARERR1 47HGENER ATE INPUT CARD HAS MORE THAN 41 ENTRIES $8 

20 K = 1 

All ) = A TE( 1) 

MM 1=M- i 

IF(MMl.EQ.O) GO TO 60 
DO 50 1=1, MM1 

SGN=SI6N< l.,ATE12*l+l)-ATE(2*I-l)) 

D=SIGNCATE<2*I),SGN) 

WRITE! 6,10031 ATE 12*1-1) , 0 , ATE ( 2 *1 +1 ) 

1003 FORMAT (45X,G13.6, 1H 1 , G 13. 6 , 1H ) , G13. 6 ) 

H=D* l . E— 4 

C=0. 

30 K =K +1 

1F1 K.GT.N )CALL ARERR1 43HGENARY IS BEING ASKED TO OVERFILL A VECTOR 

1 $) 

C=C+1. 

B=ATE( 2* I- l) +C*D 

IFlia+H )*SGN.GT.ATE(2*I+1)*SGN)G0 TO 40 
AIK )=B 
GO TO 30 

40 AIK ) = A TE( 2*1 + 1) 

50 CONTINUE 
60 L=K 
RETURN 
END 
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6 

S TO 

0 

store 


C LA 

FLC 

FLOAT 


lLS 

27 

ROUND 


F VU 

C 

H A ND U 


3 70* 

3 ? '■+ 



r ka 

1 


FLC 

OCT 

0 0 O') 000 ) 02 30 

CON S 

0 EC 

3051 75/8125 

ON fc 

3 tc 

1 


C 

OCT 

1 7 DO 00 0002 00 


END 




THE LEW ORDER PART AT DAM 
MuRMALIZE » AND 
THE 
NO. 


EXPUNENT OF RANDOM NO. 

5 EXP 15 

ONE 

NORMALIZING CONSTANT 
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